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Section  1 

OBJECTIVES 


An  analytical-numerical  study  was  pursued  by  the  present  investigators  under  AF- 
SOR  sponsorship  during  February  1990  -  May  1992.  The  primary  objectives  were: 
(i)  To  characterize  the  unsteady  separation  accompanying  dynamic  stall  of  an  airfoil 
in  constant-pitch-rate  motion;  (ii)  To  develop  a  flow-control  strategy  to  optimize  the 
instantaneous  maximum  lift  and  the  sustained  lift;  (iii)  To  develop  3-D  unsteady  flow 
simulation  capability  to  examine  flow  past  a  rectangular  planform  wing  as  well  as  a 
delta  wing.  In  connection  with  the  last  primary  objective,  the  secondary  objectives  of 
the  development  of  a  three-dimensional  unsteady  Navier-Stokes  analysis  using  itera¬ 
tive  schemes,  as  well  as  computational  efficiency  and  management  and  post-processing 
of  large  data  bases  become  very  significant. 

A  multi-tasked  research  effort  was  undertaken  to  achieve  these  objectives.  A  phys¬ 
ically  consistent  2-D  unsteady  flow  analysis  for  maneuvering  airfoils  is  developed  in 
theory  and  direct-simulation  code,  using  the  concept  of  circulation  in  a  vorticity- based 
formulation.  An  unsteady  flow-control  strategy  is  developed  for  effective  management 
of  the  associated  lift.  The  analysis  formulated  for  studying  3-D  unsteady  flow  is  also 
based  on  vorticity  and  velocity,  but  uses  an  iterative  technique.  In  addition,  to  ac¬ 
celerate  convergence,  an  efficient  multi-grid  algorithm  is  developed  for  the  velocity 
problem  emerging  from  this  formulation.  The  issue  of  efficiency,  at  least  for  2-D  flows, 
is  further  addressed  by  developing  a  temporally  adapting  grid  in  which  the  adaption 
is  based  on  the  evolving  flow  solution.  This  adaptive  grid  technique  has  the  potential 
to  resolve  very  fine  scales  that  are  critical  in  the  study  of  dynamic  stall. 

The  significant  accomplishments  made  toward  achieving  the  above-stated  objec¬ 
tives  are  briefly  described  next. 
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Section  2 


DESCRIPTION  OF 

SIGNIFICANT 

ACCOMPLISHMENTS 


All  of  the  areas  of  research  pursued  and  the  progress,  as  well  as  the  specific  achieve¬ 
ments,  made  in  these  studies  during  the  27-month  grant  period  axe  briefly  summarized 
in  the  following  sections. 

2. A  2-D  Dynamic  Stall  and  Its  Control 

2.A.1  Introduction 

The  realization  of  supermaneuverable  flight  necessitates  the  use  of  unsteady  non¬ 
equilibrium  flow  analyses  and  examination  of  a  more  comprenensive  parameter  en¬ 
velope  which  can  capture  significant  time-dependent  effects.  The  numerous  complex 
flow  phenomena  and  interactions  which  might  occur  during  a  supermaneuver  in  the 
high-alpha  (generally  between  30°  and  90°)  flight  regime  are  highly  nonlinear  pri¬ 
marily  due  to  flow  separation,  presence  of  vortex  dominated  flows  and  unsteadiness. 
As  a  consequence,  there  is  a  strong  coupling  between  the  prevailing  aerodynamics 
and  flight  dynamics,  sometimes  leading  to  chaotic  flow  and  loss  of  control  of  the  air¬ 
craft  due  to  tail  spin.  To  improve  the  maneuverability  of  high-performance  military 
aircraft  at  very  high  angle  of  attack,  researchers  are  currently  developing  analytical 
and  experimental  design  tools  that  can  take  full  advantage  of  unsteady  aerodynam¬ 
ics.  When  achieved  by  operating  in  the  post-stall  flight  regime,  supermaneuverability 
permits  improved  combat  capability. 

The  initiation  of  a  high-alpha  maneuver  may  involve  large-amplitude  deflections 
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of  control  surfaces  or  rapid  pitch-up  motions  of  the  lifting  surface  itself.  During  this 
post-stall  maneuver,  instead  of  experiencing  massive  flow  separation  and  the  result¬ 
ing  loss  of  lift,  the  lifting  surface  develops  an  energetic  dynamic-stall  vortex,  which 
temporarily  leads  to  a  significant  increase  in  lift  and  drag  forces.  Thus,  during  the 
time  that  the  dynamic-stall  vortex  is  created  on  the  suction  surface  and  convects 
over  it,  the  flow  field  is  characterized  by  a  number  of  dominant  Bow  features  which 
include  growth  of  boundary  layer,  separation,  unsteadiness,  primary  and  secondary 
shear  layer  instability,  unsteady  separation,  shock/boundary-layer  and  inviscid-viscid 
interactions,  and  vortex-vortex  and  vortex-surface  interactions.  Thus,  the  dynamic 
stall  event  is  richly  endowed  with  many  basic  fluid  mechanics  phenomena.  The  un¬ 
derstanding  and  control  of  this  event  is  important  not  only  for  supermaneuverable 
aircrafts,  but  also  for  helicopter  rotor  blades,  compressor  blades,  wind  turbines,  etc. 

Towards  the  development  of  a  supermaneuverable  flight  capability,  Lang  and  Fran¬ 
cis  (1985)  had  articulated  many  of  the  research  problems  that  are  currently  being 
pursued  by  researchers.  Many  of  the  cor  epts  discussed  by  them  are  now  being 
tested  in  the  first  international  X-plane,  namely,  the  X-31,  known  officially  as  the 
Enhanced  Flight  Maneuverability  Demonstrator;  see  Lerner  (1991).  This  is  the  first 
aircraft  which  has  successfully  maneuvered  in  the  post-stall  regime,  thereby  making 
this  liability  an  asset.  Thus,  any  insight  gained  in  understanding  the  dynamics  of 
the  post-stall  regime  could  be  very  valuable  in  improving  the  performance  and  flight 
envelope  of  this  type  of  aircraft. 

Since  April  1990,  there  have  been  three  major  workshops/conferences,  where  the 
research  issues  as  well  as  the  results  of  the  studies  in  the  area  of  post-stall  maneuvers 
have  been  discussed.  References  [A. 2,  A. 6  and  A.  16]  provide  sources  where  current 
work  is  published  and  the  reader  would  benefit  immensely  from  them.  In  light  of 
these  references,  together  with  recent  excellent  review  articles  by  Carr  (1985),  Helin 
(1989)  and  Visbal  (1990),  it  is  decided  to  not  include  a  review  of  the  literature  in  the 
present  paper. 

The  present  authors  have  been  studying  forced  unsteady  separated  flows,  using 
the  NACA  0015  airfoil  for  some  time.  Specifically,  K.  Ghia,  Osswald  and  U.  Ghia 
(1990)  provided,  for  very  low  Re,  results  that  very  vividly  showed  four  stages,  as 
classified  by  Walker  (1992),  that  describe  the  dynamic  stall  event.  These  are: 
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i.  triggering  phase,  in  which  flow  separation  occurs  near  the  leading  edge  (LE) 
on  the  suction  surface  due  to  adverse  pressure  gradient; 

ii.  separation,  in  which  vorticity  accumulates  in  the  surface  layer  and  the  onset 
of  interaction  with  the  outer  fluid  is  about  to  occur; 

iii.  strong  interaction,  in  which  a  vorticity  plane  erupts;  and  finally, 

iv.  inviscid  interaction  phase,  in  which  the  dynamic-stall  vortex  is  formed  due 
to  roll-up  of  the  free  shear  layer,  which  entrains  fluid  from  the  boundary  layer. 

The  details  of  the  analysis  were  given  by  Osswald,  K.  Ghia  and  U.  Ghia  (1990). 
Subsequently,  K.  Ghia,  Yang,  Osswald  and  U.  Ghia  (1991)  provided  the  results  for 
higher-Re  flows  by  treating  the  nonlinear  convection  terms  using  a  third-order  accu¬ 
rate  biased  upwind  differencing  scheme,  while  still  retaining  the  central  differencing 
scheme  for  all  other  spatial  derivatives.  Thereafter,  K.  Ghia,  Yang,  Osswald  and 
U.  Ghia  (1992)  demonstrated  the  suppression  of  the  dynamic-stall  vortex  using  an 
active  control  strategy  of  suction/injection.  In  addition,  K.  Ghia,  Yang,  U.  Ghia  and 
Osswald  (1992)  successfully  studied  the  dynamic-stall  phenomenon  using  a  modified 
NACA  0012  airfoil  undergoing  sinusoidal  oscillations.  The  primary  objective  in  that 
investigation  was  to  simulate  and  analyze  the  Grand  Challenge  Problem  posed  by 
Carr  (1990).  The  experimental  data  for  this  problem  was  that  of  McAllister  and 
Carr  (1979).  Further,  Osswald,  K.  Ghia  and  U.  Ghia  (1992)  extended  the  analysis  to 
provide  a  more  accurate  far-field  boundary  condition,  which  required  the  coupling  of 
viscous  circulation  T(t)  with  the  original  analysis  of  the  authors,  Ref.  [A. 13],  which 
used  thereby  arriving  at  a  (cJ/,  T(t))  formulation.  In  addition,  some 

results  were  provided  for  Re  =  45,000  with  constant  pitch-rate  dt+  =  0.2,  and  Re  = 
52000  with  a+  =  0.072. 

The  primary  objective  of  the  present  study  is  to  extend  the  analysis  of  K.  Ghia, 
Yang,  Osswald  and  U.  Ghia  (1992)  to  provide  an  accurate  far-field  boundary  condition 
by  correctly  implementing  the  prevailing  viscous  circulation  there.  Although  this 
analysis  parallels  that  of  Osswald,  K.  Ghia  and  U.  Ghia  (1992),  it  has  some  significant 
differences  in  the  way  me  mathematical  formulation  is  set  up  and  particularly  in  the 
details  of  the  numerical  procedure.  In  addition,  it  is  also  the  goal  of  this  study  to 
analyze  the  simulation  results  more  fully  in  light  of  the  available  experimental  data 
or  numerical  results,  so  as  to  assess  the  overall  2-D  (cJj,  ifrf,  T(t))  formulation,  where 
T(t)  is  the  viscous  circulation  at  infinity. 
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2. A. 2  Mathematical  Formulation 


The  time-dependent  flow  around  a  maneuvering  airfoil  is  governed  by  the  unsteady 
Navier-Stokes  (NS)  equations.  The  selection  of  the  specific  forms  of  these  nonlinear 
coupled,  partial  differential  equations  used  in  the  present  work  is  based  primarily  on 
two  factors,  namely,  (i)  dependent  variables,  and  (ii)  reference  frame. 

On  the  Choice  of  Dependent  Variables 

K.  Ghia,  Hankey  and  Hodge  (1977)  used  a  regular  grid  and  solved  the  primitive- 
variable  form  of  the  NS  equations  in  which  the  continuity  equation  was  satisfied  using 
the  Poisson  equation  for  pressure.  Osswald  (1981)  examined  various  2-D  formulations 
of  the  NS  equations  and  concluded  that  the  2-D  primitive-variable  formulation,  with 
Poisson  equation  for  pressure,  should  be  discretized  using  a  staggered  grid,  rather 
than  a  regular  grid,  if  the  discrete  equations  are  to  be  exactly  consistent  with  the 
continuous  equations.  U.  Ghia  and  K.  Ghia  (1987)  had  reviewed  various  NS  analyses 
for  3-D  flows;  their  study  included: 

i.  velocity- pressure  (V,p), 

ii.  velocity- vorticity  (V,u>)  and 

iii.  vector-potential- vorticity  (A,u>)  formulations. 

Osswald,  K.  Ghia  and  U.  Ghia  (1987)  clarified  further  that,  although  primitive- 
variable  formulations  are  widely  used  due  to  their  popularity  in  compressible-flow 
analyses,  the  velocity-vorticity  formulation  is  equally  competitive  and  perhaps  more 
advantageous  because  it  directly  provides  the  vorticity  which  is  the  most  relevant 
quantity  in  the  flow.  In  addition,  it  was  pointed  out  that  the  [V  ,u)  formulation 
leads  to  a  natural  decoupling  of  the  governing  equations,  since  the  spin  dynamics  of 
a  fluid  particle  governed  by  the  vorticity-transport  problem  can  be  decoupled  from 
the  translational  kinematics  of  the  fluid  particle  represented  by  the  elliptic  velocity 
problem.  This  is  not  the  case  for  the  primitive-variable  formulation.  It  was  also 
pointed  out  that  a  careful  analysis  of  the  (V,u>)  formulation  could  reduce  its  com¬ 
putational  requirements  to  match  those  of  the  (V,p)  formulation.  Finally,  it  was 
pointed  out  that  the  vector-potential-vorticity  formulation  required  specification  of  a 
non-physical  boundary  condition  to  correctly  set  up  the  problem  and  that  this  poses 
numerical  difficulties.  Huang,  U.  Ghia  and  K.  Ghia  (1992)  have  revisited  this  issue 
and  have  provided  some  additional  details.  Gatski  (1991)  as  well  as  Gresho  (1991) 


5 


have  also  reviewed  the  various  NS  formulations  and  have  provided  some  insight  into 
the  selection  process  for  the  dependent  variables. 

On  the  selection  of  Reference  Frame 

Even  for  maneuvering  bodies,  it  is  possible  to  work  with  an  inertial  reference  frame, 
with  boundary- aligned  coordinates  being  computed  at  every  instant  of  time  to  provide 
for  the  body  motion.  Chyu  (1981)  as  well  as  Saiari  and  Roache  (1990)  have  provided 
analyses  using  an  inertial  reference  frame.  The  major  advantage  with  this  approach 
is  that  the  far-field  boundary  conditions  remain  undisturbed.  On  the  other  hand, 
Mehta  (1977),  Sankar  and  Tassa  (1980),  Visbal  and  Shang  (1989)  and  the  present 
authors  have  elected  to  work  with  the  body-fixed  non-inertial  reference  frame.  In  this 
approach,  although  the  grid  remains  undistorted,  the  far-field  boundary  conditions 
require  special  attention.  In  the  present  study  of  2-D  forced  unsteady  separated  flow, 
the  ( V,u> )  formulation  is  used  in  generalized  coordinates  in  the  body- fixed  non-inertial 
reference  frame.  Osswald,  K.  Ghia  and  U.  Ghia  (1990)  have  shown  that,  for  the 
( V,Q )  formulation  with  divergence  and  curl  operators  expressed  in  the  generalized- 
coordinate  non-inertial  reference  frame,  the  inertial  vorticity  diffuses  as  for  the  case 
of  a  fixed  body,  but  advects  with  apparent  velocity  rather  than  with  inertial  velocity. 
Thus,  except  for  the  appearance  of  the  apparent  velocity,  the  velocity-vortidty  for¬ 
mulation  of  the  unsteady  NS  equations  is  “nearly”  form-invariant  under  a  generalized 
non-inertial  coordinate  transformation.  This  offers  a  significant  advantage,  in  that  it 
leads  to  a  unified  algorithm  for  both  non- maneuvering  and  maneuvering  body  flows. 
Thus,  in  the  present  study  for  2-D  flow  past  a  maneuvering  body,  the  (i/>,  Q)  formula¬ 
tion  is  used  in  a  body-fixed  non-inertial  reference  frame.  This  form  has  been  shown 
by  Osswald,  K.  Ghia,  and  U.  Ghia  (1988)  to  be  equivalent  to  the  {V,u)  formulation, 
but  it  is  computationally  more  efficient.  Speziale  (1987)  had  also  shown  that,  for 
the  special  case  of  rotation,  the  (V^w)  formulation  is  form- invariant  and  that,  under 
this  condition,  non-inertial  effects  will  enter  the  solution  only  through  the  initial  and 
boundary  conditions. 

Governing  Differential  Equations 

An  arbitrary  maneuver  can  be  completely  defined  by  specifying  the  trajectory  rs//(0 
of  some  point  B  fixed  on  the  body,  (see  Fig.  A.l),  with  respect  to  an  inertial  observer, 
together  with  the  specification  of  the  instantaneous  angular  velocity  Qg(f)  of  the  body. 
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Kinematically,  the  translational  velocity  of  the  origin  B  of  the  body-fixed  frame  is  then 
Ve//(t)  —  d(rB/i)/dt,  and  the  translational  acceleration  is  dg/iit)  =  iP(fg/j)fdt 2, 
while  the  angular  acceleration  is  ds(t)  =  d(ftB)/dt.  In  the  present  analysis,  these 
functions,  which  define  a  specific  maneuver,  are  assumed  to  be  explicitly  prescribed 
functions  of  time,  and  will  be  given  in  a  later  sub-section. 

For  an  inertial  reference  frame,  the  unsteady  incompressible  Navier-Stokes  equa¬ 
tions  are  given  in  terms  of  the  derived  variables  (  V/,u7/  )  as 


Continuity  Equation  and  Kinematic  Definition  of  Vorticity  : 


V  /  •  V/  =  0  , 

(A.l) 

V  i  x  Vj  —  <jj[ . 

(A.2) 

Vorticity- Transport  Equation  : 

~  +  V,  x  (0,  x  V,)  +  i(V,  x  V,  X  a,)  =  0 , 

(A.3) 

where  Re  —  such  that  is  the  reference  free-stream  velocity  and  c  is  the  airfoil 
chord  length.  Also,  the  subscript  I  on  V/  denotes  that  the  implied  spatial  differen¬ 
tiation  is  with  respect  to  inertial  coordinates.  The  transformation  of  Eqs.  (A.1-A.3) 
to  a  generalized  coordinate  non-inertial  body-fixed  reference  frame,  as  carried  out  by 
Osswald,  K.  Ghia  and  U.  Ghia  (1990),  leads  to  the  following  form: 


Continuity  Equation  and  Kinematic  Definition  of  Vorticity  : 


<1 

II 

o 

(A.4) 

VxV/  =  u/. 

(A.5) 

Vorticity- Transport  Equation  : 

Ii  +  Vx(i,xi')+  j--(V  x  V  x  a,)  =  0, 

(A.6) 

where  all  divergence  and  curl  operators  are  with  respect  to  the  generalized-coordinate 
non-inertial  reference  frame,  the  apparent  velocity  V  is  kinematically  related  to  the 
inertial  velocity  as 


V  =  Vi  -  Vfl//(f)  -  ClB(t)  x  f,  (A. 7) 
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and  the  position  vector  f  is  given  as 


f  -  rt  ~  rBfi(t) ,  (A.8) 

with  Vs//(i),£ij3(f)?fB//( t)  being  the  known  functions  which  define  a  specific  ma¬ 
neuver. 


Hence,  except  for  the  appearance  of  the  apparent  velocity,  the  velocity-vorticity 
formulation  of  the  unsteady  Navier-Stokes  equations  is  form-invariant  under  a  gen¬ 
eralized  non- inertial  coordinate  formulation.  Therefore,  the  algorithm  developed  to 
solve  Eqs.  (A.1-A.3)  remains  valid  for  the  solution  of  Eqs.  (A.4-A.6).  The  analysis 
described  so  far  is  valid  for  3-D  flows  with  six-degree-of- freedom  maneuvers. 


For  the  2-D  simulations  of  interest  here,  it  is  computationally  more  efficient  to  use 
the  variables  rather  than  the  (V\to)  variables.  Further,  due  to  the  unbounded 

nature  of  the  stream  function  in  the  far  field  at  true  infinity,  it  is  essential  that  a 
deviational  stream  function  be  employed  as  the  dependent  variable;  this  is  given  as 


0/  =  y  ~  VB/i(t )  +  rp? t)  +  ^  Q  (A-9) 

where  y  —  yB/i(t)  is  the  vertical  inertial  coordinate  passing  through  the  instantaneous 
location  of  the  pitching  axis.  The  viscous  circulation  T(t)  is  positive  in  the  clockwise 
direction  and  ‘a’  is  the  radius  of  the  circle  to  which  the  airfoil  is  transformed.  The 


inertial  velocity  then  becomes 


V{  =  i  + 


Real(Vmei) 

9n 


+ 


i  w 

V9dt2. 


ei  + 


RealjVh) 

912 


i  w 


e2. 


(A. 10) 


Here,  V’  is  the  complex  velocity  defined  as  i~^  due  to  the  generation  of  circulation, 
and  i  is  the  unit  base  vector  of  the  inertial  Cartesian  reference  frame  shown  in  Fig. 
A.l.  Further,  ei  and  e2  are  the  covariant  base  vectors  of  the  generalized-coordinate 
(£*,  £2)  non-inertial  body-fixed  reference  frame.  The  determinant  of  the  metric  tensor 
g  appearing  in  Eq.  (A.  10)  is  given  as 


9  —  {9u9n  ~  9u9ii) 


where 


The  independent  variables 
whose  unit  base  vectors  are 


2  (dxk\  (dxk\ 

9“~h  Uv  U'/ 


(A. 11) 


(xl,  x 2)  represent  the  body-fixed  Cartesian  coordinates 
and  e2  as  shown  in  Fig.  A.l. 
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In  light  of  Eqs.(A.7)  and  (A. 10),  the  contravariant  components  of  the  apparent  ve¬ 
locity  are 


{sin  0B(t)  -  Vg/t(t)  -  x lnB(t)}e2  ■  el  +  Real{Ve,)}  +  (A.12) 

and 

W  =  ^[{cos ^s(t)  —  V^/;(t)  4- i2flB(t)}e i  ■  e2 + 

fl.l.D 

(sin^fl(t)  —  Vjyj(t)  —  ‘  ^2  +  fieal(V*e2)]  — ,(A.13) 

where  (-0s(t)]=  af(t)  is  the  instantaneous  pitch  angle  of  the  airfoil  and 
V B/l(t )  =  i  +  ^b//(0c  2- 


Navier-Stokes  equations  in  generalized  orthogonal  non-inertial 
body-fixed  coord  mat esytake  the  following  for 
Stream  Function  Equation: 


Vorticity  Transport  Equation: 

W  +  Mv^l)  +  ^2  («H^)  = 

1_  /  5  /y22  dwA  ] 

/?e\^  ^aev  **  wS^vJ* 


(A.14) 


(A. 15) 


In  Eqs.  (A. 14)  and  (A. 15),  the  elliptic  problem  for  the  deviational  stream  function  is 
coupled  with  the  temporally  parabolic  vorticity- transport  equation.  The  mathemati¬ 
cal  formulation  can  be  completed  by  providing  the  boundary  and  intitial  conditions. 
However,  before  discussing  these,  the  control  strategy  developed  for  this  study  will  be 
examined.  The  boundary  and  initial  conditions  can  then  be  presented  in  their  final 
form. 


Active  Control  Using  Suction  and  Injection 

From  a  theoretical  consideration,  Osswald  (1992)  has  pointed  out  that  the  vorticity 
created  at  the  surface  drives  some  of  the  unsteady  flow  phenomena.  Thus,  flow  phe¬ 
nomena  such  as  incipient  separation  and  subsequent  development  of  vortical  structure 
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near  the  surface  perhaps  can  be  controlled  by  managing  the  creation  of  surface  vortic- 
ity,  which  is  governed  by  the  vorticity  boundary  condition  at  the  surface  and,  hence,  is 
external  to  the  governing  equation.  Thus,  modulated  suction/injection(MSI)  control 
is  developed.  The  active  control  model  developed  is  based  on  the  principal  objective 
of  delaying  the  unsteady  separation  to  suppress  the  formation  of  the  dynamic-stall 
vortex.  A  secondary  goal  is  to  increase  the  lift  and  reduce  drag.  In  additition  to  these 
objectives,  the  constraints  used  in  developing  the  active  control  model  are  that 

i.  the  maximum  non-dimensionalized  main  flow  rate  be  atmost  1%, 

ii.  the  net  mass  addition  be  zero  and,  finally, 

iii.  the  control  model  for  MSI  be  as  simple  as  possible. 


Constraint  (ii)  is  not  critical;  however,  if  employed,  it  allows  the  downstream 
boundary  condition  for  the  basic  flow  without  MSI  to  remain  valid  for  the  flow  even 
with  the  implementation  of  active  control. 

The  active  control  model  depends  on  a  large  number  of  variables,  some  of  which 
are  location,  magnitude,  suction/injection  rate,  duration,  etc.  A  careful  numerical 
optimization  led  to  a  trapezoidal  profile  for  the  MSI  velocity,  with  the  key  parameters 
being  depicted  in  Fig.  A. 22.  In  the  present  study,  the  suction  velocity  is  taken  to 
be  between  3.5%  and  4.5%  of  the  free-stream  velocity,  and  is  applied  on  the  upper 
surface  of  the  airfoil  over  approximately  9%  of  the  chord  starting  from  the  leading  edge 
(LE).  The  physical  injection  velocity  is  slightly  smaller  in  magnitude  as  compared  to 
the  physical  suction  velocity.  It  should  also  be  pointed  out  that  the  segment  of  the 
airfoil  surface  in  the  computational  plane  over  which  suction  is  applied  is  the  same  as 
that  over  which  injection  is  applied.  Then,  to  satisfy  the  constraint  of  zero  net  mass 
addition,  injection  is  applied  over  14%  of  the  chord  on  the  lower  surface  near  the 
trailing  edge  (TE)  of  the  airfoil.  In  fact,  the  constraint  can  be  expressed  as  follows. 
In  the  physical  plane,  for  suction  applied  over  the  upper  surface  near  the  LE  segment 
(a-b),  and  injection  over  the  lower  surface  near  the  TE  segment  (c-d),  zero  net  mass 
addition  requires 

I"  ~(V  ■  it)il  U,  +  ld  ~(V  ■  e3)dl  U  =  0  ,  (A. 16) 

Ja  \J922  \J  922 

where  the  apparent  velocity  V  on  the  surface  can  be  written  as 

V  =  =  T({1, £,!)-§=  +  N((\ .  (A.17) 

y/9u  \Zg-i2 
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In  the  non-inertial  reference  frame  of  the  generalized-coordinates  (£\£2),  Eq.  (A. 16) 
is  given  as 

f  y/^me,(it)de=o.  (a.is) 

np 

where  represents  the  non-porous  segments  of  the  airfoil  surface  in  the  computa¬ 
tional  plane. 


Boundary  Conditions 

Using  Eqs.  (A. 12)  and  (A. 13),  the  apparent  velocity  V  can  be  expressed  as 

V  =  jcos0a(O  —  Us//(t)  —  x2flfl(t)j  +  |sin(?a(t)  —  Uj//(t)  4-  xlOB(t)]  e 2 

,  RealiV'er)  ,  1  .  ,  f/?ea/(U*e2)  1  d0fl  .  _ 

+  - +  ~P~Z77  ei  + - e2- 

£11  v/5  l  g22  V?  df1 

From  Eqs. (A.  17)  and  (A. 19), 

l6o^=  Kcos  MO  -  Vb/M  +  *J*M0)«i  + 

{sin0B(O  -  Uj/jr(f)  -  xlftB(t)}e2]  •  e2 

+i?eai(V*e-2)  -  v®V(£l,  t) .  ( A.20) 

Eq.  (A.20)  leads  to  the  boundary  condition  for  the  stream  function  at  the  surface  as 

(fSfi.O  =  {*TC//(0  “  «*M01  -  sl[VB//(0  -  sin  MO] 

-  j x  ,  (A.21) 


where  <f*p  represents  the  non-porous  segments  of  the  body  surface.  Now,  the  surface 
vorticity  can  be  evaluated  from  the  following  equation: 


.  _ 1  \d  (gnd^?\  9  /gug^fV 

vs  kl  9t‘ )  se  \^g  de  ’ 

subject  to  the  constraint  of  zero  slip,  which  leads  to 

Iw,=  -[{C0»#B(()  -  t)  +  rJ!lB(0}«i  + 
{sin0B(t)  -  V^jit)  -  x^aCO}^]  •  et 
—Real(V*ei)  +  >/5ur(£l»£>0 . 


(A.22) 


(A. 23) 


Eq.  (A.20)  is  also  used  in  the  evaluation  of  Eq.  (A.22).  The  corresponding  condition 
in  the  far  field  can  be  written  as 

=  0,  (A. 24) 


n 


Ut  =  0.  (A. 25) 

The  examination  of  this  boundary  condition  reveals  that  the  viscous  circulation  T(t) 
is  still  not  determined;  its  determination  is  described  in  the  next  sub-section. 


Viscous  Circulation  and  Necessary  Condition  for  Its  Determination 


The  inviscid  circulation  r,n(t)  is  given  as 


rin(<)  =4*-a  sin  a/(t)  (A. 26) 

where  ‘a’  is  the  radius  of  circle  to  which  the  airfoil  is  transformed(see  Fig.  A.3(d,  e)) 
and  q/  is  the  flow  angle  at  the  given  instant  of  time.  The  circulation  T(t)  for  the 
overall  viscous  flow  is  related  to  the  inertial  velocity  V}  as 

r(0  =  /  Vfdl=-  Vr  dl,  (A. 27) 

J  X>  JC00 

taken  around  a  circle  C '<*>  with  its  radius  equal  to  infinity.  The  use  of  the  Stokes’ 
theorem  leads  to 


r(t)  =  —  /  (V  x  Vi)  ■  nds  =  —  /  £/  •  nds, 
Js0  Js0 


(A.28) 


where  Sa  is  the  corresponding  area  of  the  entire  flow  domain  whose  boundary  curve 
is  Coo.  Even  with  this  condition,  the  stream- function  and  vortidty  equations,  Eqs. 
(A.  14- A.  15),  still  constitute  a  singular  system  of  equations  due  to  the  linear  depen¬ 
dence  amongst  these  equations.  This  led  to  imposing  the  condition  that  the  pressure 
in  the  multiply-connected  flow  domain  be  single  valued  and  continuous  along  the 
body  surface,  i.e., 


dl  =  0. 


(A. 29) 


For  flow  past  a  square  cylinder  in  a  channel,  again  a  multiply-connected  domain, 
Matlda,  Kuwahara  and  Takami  (1975)  had  used  Eq.  (A. 29)  to  close  the  equation  set 
for  their  internal  flow  problem.  This  has  provided  the  impetus  to  use  this  equation 
in  the  present  analysis  also.  Indeed,  Eq.  (A. 29)  is  used  in  the  present  analysis 
and  serves  to  close  the  equation  set.  This  condition  is  an  alternative  statement  of 
the  Kutta  condition  for  viscous  unsteady  flow,  namely,  that  the  pressure  at  the  TE 
has  to  be  continuous  and  single  valued.  Osswald,  K.  Ghia  and  U.  Ghia  (1992)  have 
pointed  out  that  the  only  influence  that  pressure  has  on  an  incompressible  flow  is  that 
it  controls  the  creation  and  distribution  of  surface  vorticity.  Further,  this  condition 
helps  to  close  the  problem  thereby  implying  that  advanced  wall  vorticity  can  now  be 
accurately  determined.  Next,  the  initial  conditions  used  in  the  present  analysis  are 
discussed. 
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Initial  Conditions 


To  determine  the  flow  field  past  a  maneuvering  body,  initially  the  simulation  of  the 
flow  past  a  stationary  airfoil  at  fixed  angle  of  attack  is  carried  out.  For  this  latter 
stationary  configuration,  the  motion  is  started  impulsively  from  rest  at  t=0  and  he 
corresponding  inviscid  solution  is  used  as  the  initial  condition.  The  viscous  flow  past 
a  stationary  body  at  a  given  angle  of  attack  is  determined.  Next,  up  to  an  asymptotic 
time,  f  =  T,  the  maneuvering  motion  is  initiated  at  i  =  T  and  the  viscous  circulation 
at  the  start  is  taken  as 

r(«)  u=  r0(r)  (a. 30) 

where  limt-_2,r(t)  =  ro(T)  is  the  asymptotic  value  of  the  circulation  computed  by 
accounting  for  the  angular  acceleration  Qs//(t)  due  to  manuevering  motion  in  the 
value  of  the  circulation  around  a  stationary  airfoil. 


Pressure  and  Force  Coefficients 

The  analysis  to  be  discussed  here  follows  the  development  provided  by  Yang  (1992). 
In  the  present  (u>/,  ,  T(t))  formulation,  pressure  does  not  appear  and,  as  such,  to 

determine  pressure,  it  is  necessary  to  revert  back  to  the  linear-momentum  equation 
given  as 

|  +  (Vxl/,)xi';  +  i(VxVxV))s-v|p+  .  (A. 31) 

The  pressure  field  in  the  incompressible  flow  can  be  evaluated  accurately  using  the 
pressure  Poisson  equation  obtained  by  taking  the  divergence  of  the  linear-momentum 
equation,  Eq.  (A. 31).  If,  on  the  other  hand,  only  surface  pressure  is  desired,  it  can  be 
determined  more  readily  by  carefully  evaluating  Eq.  (A. 31),  along  the  surface;  this 
leads  to 


dip  Hr  <?)  , 

I  body 


r-dT(e,C,t)  -x  - 

-y/9n - ^ -  -{aB/[[t)  X  n  ej 

s - v— — - - ' 

Tangential  Acceleration 
due  to  circulation  control 
along  the  body  surface 

Pressure  gradient  due  to 
mass  transfer  along  the 
body  surface  as  a  result 
of  MSI 


(A. 32) 
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where 

q  =  5  {TH(\ (1,0  +  ^({'.ft.O)  +“'bM‘>x'  +  4//(0*j  -  \  {(x'f  +  (4)’}  4(0- 

"  ' - - — ^ - ' 

Addition  of  kinetic  energy  due  to 
MSI 

(A. 33) 

For  evaluation  of  the  terms  involving  apparent  velocity  on  the  surface,  not  all  terms 
are  zero  and,  for  some  additional  details,  see  Yang  (1992).  Careful  examination 
of  Eq.  (A. 32)  reveals  that  the  temporal  increase  of  the  tangential  control  velocity 
can  accelerate  the  flow,  since  <  0  and,  for  suction  with  t)  <  0,  as 

implemented  in  the  unsteady  separation  region  where  uj  >  0,  y/gnNd1,^,  t)u>i  will 
also  produce  favorable  pressure  gradient  and  will  lead  to  reduction  in  the  boundary 
layer  thickness.  Thus,  integrating  Eq.  (A. 32)  leads  to 


P  l£f=  -q  If?;  +  /f  (v5^VK‘,*i,0u,/  -  V3I7 

■'si  ^ 


"ij‘  X  f) '  6‘‘  +  Te^w}^' '  <A'34) 

This  then  permits  the  determination  of  the  surface  pressure. 

The  force  and  moment  coefficients  as  well  as  the  surface  force  due  to  viscosity  are 
now  calculated  using  the  surface  pressure  from  Eq.  (A. 34),  and  the  definition  of  the 
coefficient  of  pressure  as 


c’-'ipui-2p’ 

where  p'  is  the  dimensional  pressure. 

Further,  the  force  coefficient  Cp  is  defined  as 


(A.35) 


r  x  _ 
LFP  -  1 


pULcb  ’ 


where  the  spanwise  length  b  is  taken  as  unity. 

Substituting  Eq.  (A. 34)  into  Eq.  (A.35),  Eq.  (A. 36)  is  then  written  as 


CF,  =  -  /  (  ^Cre 4  a'  . 

Jbody  ( \fgri  ) 


Similarly,  the  net  dimensional  force  Fv  due  to  viscosity  is  given  as 


(A. 36) 


(A. 37) 


K  =  j  2p  (I  [(v-q-)  +  q-v-)]  - 1 .  h<u- ,  (a.3s) 
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where  jxm  is  the  dynamic  viscosity  and  /  is  the  identity  tensor. 
The  coefficient  of  the  viscous  shear  force  is  defined  as 

n  ^ v 

=  — 


(A. 39) 


\puicd 

Use  of  the  expression  for  apparent  velocity,  Eq.  (A. 7),  permits  Eq.  (A. 39)  to  be 
expressed  as 


Cr-=kL^V)  +  (VV)  I'"*' 


(A. 40) 


For  flows  with  MSI,  an  expression  for  Eq.  (A. 40)  can  be  obtained  in  non-inertial 
body-fixed  generalized  coordinates  (£\<f2)  as 

dV  dV 


-  2  r  \^dV  _2dV 

CFu  ~  Re  Ldy  \6  d^+e  df,2  + 


\Jjri 


-e'  >  •  e2-i7 —d£  , 


(A.41) 


where  e1  and  e2  are  the  contravariant  base  vectors  of  the  generalized  coordinates 
(4U42).  In  addition,  by  defining 


dV 

di 1 


=  Alet  +  A2e2  and 


dv 

d? 


=  Ble\  +  B2e2, 


(A.42) 


Eq.  (A.41)  can  be  rewritten  as 
9 


cv„  = 


Re 


i  +  ^B1)el+2^B2e7\de  ■  (A.43) 

2 body  {  \/g II  y/g 22  y/922  ) 


Along  the  non-porus  body  surface,  Eq.  (A.41)  is  simplified  to  give 

2 


Re 


(A. 44) 


A  similar  analysis  leads  to  the  moment  coefficients  Cmp  and  C.\fv  as 


Cm  —  ^Lcv{C)  { ^ 1  ( e i  x  e2)  +  x2(e2  x  e2)j 

2 body  [v^22  ^  J 


dt\ 


(A. 45) 


and 


C\iv  — 


/ 

Jbo 


f  v</^22 /U  ^t^Bv  1  {xJ(ei  x  et)  +  x2(e2  x  ex 


Re  Jbody  [  {  s/g^l  y/g^ 
_4_ 

~  Re 


)} 


de 


£  —^Lb1  xl(e.i  x  e2)  +  x2(e2  x  e2)|  <i£l  .  (A. 46) 

Jbody  \J922  1 


For  non-porus  body  surface, 


Cm-  =  -feL,  ^‘(i'  X  ",)  +  ^  X  ^ 


iC- 


15 


(A- 47) 


To  calculate  the  lift  and  drag  coefficients  from  Cpp  and  CVv,  it  is  convenient  to  define 
the  force  vector  F  as  a  complex  variable  such  that 


F  =  Real(F)ei  -f  Imag(F)e2  ,  (A. 48) 

With  this  definition  in  Eq.  (A. 44),  as  shown  in  Fig.  A. 2,  the  lift  and  drag  coefficients 
become 

Clp  =  Imag(CFp )  x  cos0s(t)  -  Renl(Cfp)  x  sin0a(*)> 

CDp  =  Real{CFp)  x  cos  0g(t)  +  / mag(Cpp )  x  sin^s(O) 

Ci „  =  Imag(CF„ )  x  cosOB{t)  -  /?ea/(C>v)  x  sin0a(t), 

=  Real(CFv)  x  cos  +  /ma5(Cir  )  x  sin0a(f).  (A. 49) 

Further, 

Cl  =  Clp  +  Cl ¥ , 

Cd  =  Co,  +  Cn, , 

Cat  =  C.Up  +  Cm,  • 

2.A.3  Analytical  Grid  Generation  Technique 

The  mathematical  problem  formulated  in  Section  A. 2  is  governed  by  a  set  of  partial 
differential  equations  which  can  be  solved  efficiently  using  an  advanced  numerical 
technique.  The  overall  solution  technique  will  require  a  carefully  generated  grid  that 
satisfies  numerous  criteria  as  far  as  quality  of  the  grid  is  concerned.  The  grid  should 
also  provide  the  resolution  of  dominant  length  scales  in  the  various  critical  regions 
that  evolve  in  the  flow  field.  The  use  of  appropriate  grid  topology  with  a  given 
mathematical  formulation  is  also  critical,  since  certain  Navier-Stokes  operators  may 
turn  out  to  be  singular  on  a  specific  grid  topology.  The  generalized  Schwarz- Christoffel 
grid  generation  technique  originally  described  by  Davis  (1979)  and  further  developed 
by  Osswald,  K.  Ghia  and  U.  Ghia  (1989)  for  flow  past  an  arbitrary  airfoil  is  used  in 
this  study. 

Conformal  Transformation 

The  generalized  mapping  technique  used  here  is  conformal  and  provides  analytic  met¬ 
ric  information,  guarantees  orthogonality,  and  has  the  capability  of  handling  arbitrary 


(A. 50) 
(A- 51) 
(A. 52) 
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bodies  with  sharp  corners  that  occur  for  example  at  TE,  etc.  As  shown  in  Fig.  A. 3(a), 
a  conformal  mapping  from  the  physical  z-plane  permits  the  airfoil  to  be  transformed 
into  a  line  segment  in  the  upper  half  (-plane,  as  shown  in  Fig.  A. 3(b).  A  second 
conformal  transformation  permits  this  line  segment  in  the  upper  half  (-plane  to  be 
transformed  to  the  complex  potential  plane  P,  as  shown  in  Fig.  A. 3(c).  This  is  sub¬ 
sequently  transformed  using  the  Joukowski  transformation  to  a  circle  in  the  Z-plane, 
as  shown  in  Fig.  A. 3(d).  For  lifting  cases  at  non-zero  incidence  cr/,  the  Z-plane  can 
be  rotated  by  the  effective  angle  of  attack,  ae  —  atf  +  (3  to  the  circle- R  plane,  as 
shown  in  Fig.  A. 3(e),  where  ‘a’  is  the  flow  incidence  angle  and  -/?  is  the  angle  of  zero 
lift.  The  complex  potential  P  for  a  uniform  stream  at  incidence  aj  past  an  arbitrary 
airfoil  is  obtained  from  the  knowledge  of  the  complex  poteitial  for  flow  past  a  circular 
cylinder.  The  complex  potential  P-plane  is  depicted  in  Fig.  A. 3(f)  and,  if  used  as  a 
final  plane,  it  leads  to  a  mesh  with  H-grid  topology  in  the  physical  plane.  This  H- 
type  grid  does  not  provide  sufficient  resolution  in  the  rounded-LE  region.  Therefore, 
yet  one  more  parabolic  conformal  transformation  is  used  to  transform  the  airfoil  to 
the  complex  //-plane,  as  shown  in  Fig.  A. 3(g).  This  now  leads  to  a  C-grid  topology 
in  the  physical  plane.  Even  though  the  C-grid  will  provide  some  clustering  the  LE 
region,  the  far-field  boundary  condition  to  be  placed  at  infinity  will  lead  to  numerical 
difficulties  if  7/-plane  is  used.  Hence,  the  final  transformation  contracts  this  //-plane 
to  the  computational  plane  and  is  described  in  the  next  sub-section. 

Clustering  Transformations 

One-dimensional  clustering  transformations  are  used  to  transform  the  upper-half  stag¬ 
nation  point  //-plane  tc  a  unit  square  in  the  computational  (-plane.  Osswald,  K.  Ghia 
and  U.  Ghia  (1989)  as  well  as  Rohling  (1991)  have  provided  the  details  of  the  general 
cubic-spline  functions  used.  In  the  streamwise  direction,  independent  cubic-spline 
functions  are  used  to  achieve  streamwise  grid  clustering  along  the  suction  and  pres¬ 
sure  surfaces  and  in  the  mid-wake  and  near-wake  regions,  both  above  and  below  the 
streamline  which  emanates  from  the  trailing  edge.  The  far-wake  regions  above  and 
below  the  TE  streamline  are  contracted  using  1-D  inverse  tangent  transformations. 
Similarly,  1-D  cubic-spline  transformations  are  also  used  in  the  normal  directions  in 
two  zones,  namely,  (i)  boundary  layer  and  (ii)  massively  separated  region  which  also 
includes  the  LE  shear  layer.  Finally,  a  1-D  inverse  tangent  transformation  is  again 
required  to  map  the  far-field  infinity  boundary  to  a  unit  distance  in  the  computa¬ 
tional  (-plane.  The  various  regions  discussed  here  are  depicted  in  Fig.  A. 4.  A  typical 
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grid  distribution,  for  a  NACA  0015  airfoil  with  (444  x  101)  mesh  points  is  shown 
in  Fig.  A. 5.  An  attempt  has  been  made  to  resolve  as  many  of  the  dominant  scales 
as  possible;  these  include  the  boundary  layer  scale  of  Q(Re~1^2)  for  attached  flow, 
the  streamwise  and  normal  scales  of  0(Re~ 3^8)  and  0(/?e~5/'8),  respectively,  near  the 
separation  point,  the  separated  free  shear  layer  with  the  scale  of  the  order  of  its 
thickness,  the  massively  separated  zone  of  0(1)  and  the  scale  of  0(/Ze'5/u)  for  the 
boundary  layer  eruption  in  the  separated  region.  This  is  important  in  order  to  ensure 
that  the  simulation  results  truly  preserve  the  physics  of  the  problem.  The  way  in 
which  this  is  partially  achieved  is  by  generating  a  grid  with  zonal  attributes  with 
continuous  metrics  across  all  of  the  zones  both  in  streamwise  and  normal  directions. 
Not  only  is  the  clustered  conformal  grid  generated  here  by  an  analytical-numerical 
procedure,  hut  the  evaluation  of  the  metric  tensor  and  base  vectors  is  also  carried  out 
analytically  once  the  transformations  have  been  defined. 

2. A. 4  Numerical  Method 

The  basic  numerical  technique  for  maneuvering  motion  was  developed  by  Osswald,  K. 
Ghia  and  U.  Ghia  (1990).  Variations  of  this  approach  have  been  used  by  the  present 
researchers  in  some  of  their  earlier  papers.  It  is  a  fully  implicit  finite-difference  method 
involving  direct  matrix  inversion  methodology.  All  of  the  spatial  derivatives  are  ap¬ 
proximated  using  second-order  accurate  central  difference  approximations,  except  for 
the  convective  terms.  These  latter  nonlinear  terms  are  approximated  using  a  bi¬ 
ased,  third-order  accurate,  upwind  differencing-scheme  to  be  able  to  simulate  higher 
Reynolds  number  flows.  The  vorticity  transport  equation,  Eq.  (A. 15),  is  solved  using 
a  variation  of  the  alternating-direction  implicit  (ADI)  technique  of  Douglas  and  Gunn 
(1964).  When  this  technique  is  cast  in  delta  form,  it  is  computationally  efficient.  The 
elliptic  stream  function  equation,  Eq.  (A.  14),  is  solved  by  a  direct  block  Gaussian 
elimination  (BGE)  technique.  No  explicit  artificial  dissipation  is  added  in  the  imple¬ 
mentation  of  the  overall  ADI-BGE  method. 

The  algorithm  consists  of  a  forward  elimination  sweep  for  the  vorticity  and  stream 
function  equations.  This  is  followed  by  the  determination  of  advanced  wall  vorticity, 
together  with  the  computation  of  the  viscous  circulation.  This  is  achieved  with  the 
help  of  Eq.  (A. 29)  which  requires  continuity  of  the  pressure  at  the  trailing  edge. 
Both  the  wall  vorticity  and  T(f)  are  determined  using  an  iterative  technique  and  by 
maintaining  a  strict  convergence  criterion,  such  that  the  relative  percentage  error  is 
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less  than  IQ-6,  and  the  residue  in  the  instantaneous  wall  vortidty  is  negligible.  Sub¬ 
sequently,  the  back  substitution  sweep  for  the  vortidty  and  stream  function  equations 
is  carried  out  to  achieve  the  final  solution.  It  should  be  noted  that,  due  to  the  use  of 
the  third-order  biased  upwind  differencing,  the  discretized  problem  leads  to  a  penta- 
diagonal  matrix  which  can  be  solved  using  a  generalized  Thomas  algorithm. 

One  final  comment  corresponds  to  the  determination  of  the  pressure  field.  The 
analysis  provided  so  far  is  adequate  if  only  surface  pressure  is  of  interest.  On  the 
other  hand,  if  the  entire  pressure  field  is  desired,  the  solution  of  a  Neuman:  t  Poisson 
problem  for  pressure  can  be  carried  out  as  given  by  K.  Ghia,  Yang,  Osswald  and  U. 
Ghia  (1992a).  The  results  obtained  in  the  present  study  are  discussed  next. 

2. A. 5  Results  and  Discussion 

The  continuous  and  discrete  analysis  developed  in  the  earlier  sub-sections  is  used 
to  study  the  motion  past  maneuvering  NACA  airfoils.  In  this  sub-section,  first  the 
constant-rate  pitch-up  motion  studied  will  be  outlined,  followed  by  validation  and 
verification  studies.  Finally,  some  detailed  results  will  be  provided  for  constant- 
rate  pitch-up  motion;  these  results  are  computed  without  and  with  modulated  suc¬ 
tion/injection.  Simulation  results  are  also  obtained  in  the  form  of  a  color  video  ani¬ 
mation  and,  while  discussing  the  results,  the  present  authors  have  taken  the  liberty 
to  use  these  animations  to  clarify  various  research  issues  and  provide  further  insight 
into  the  prevailing  flow  phenomena. 

A  total  of  three  flow  configurations  have  been  studied; 
these  are  as  listed  below. 


Configuration 

Re 

a+/K 

Remarks 

I 

5,000 

1.0 

for  verification  study 

II 

45,000 

0.2 

without  flow  control 

III 

45,000 

0.2 
■  .....  ■ 

with  modulated  suction/injection 

In  addition,  results  have  been  obtained  for  a  flow  configuration  with  Re=52,000 
and  a+  =  0.072.  These  results  have  been  discussed  by  Yang  (1992). 
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Constant-Rate  Pitch-Up  Manuever 

In  all  the  earlier  studies  by  the  authors,  the  constant-rate  pitch-up  motion  used  is 
the  one  defined  by  Visbal  and  Shang  (1989)  as: 
n?//(0  =  0, 


Ml)  =  ».1‘  -  (^y(l  -  (A.53) 

where  =  a+  =  ^  is  the  non-dimensional  pitch  rate  and  #e(*)  i3  the  instantaneous 
pitch  angle  of  the  airfoil.  In  the  present  study,  t0  =  0.5,  Q0  =  0.2  and  the  airfoil  is 
pitched  about  the  axis  through  the  quarter-chord  point  (QCP). 

Verification  Study 

The  earlier  studies  of  the  authors  in  Refs.  [A. 11,  A. 12,  A. 13,  A. 14  and  A. 28]  have  al¬ 
ready  established  that  the  mathematical  model  is  accurate  and  that  the  modification 
in  the  far-field  boundary  condition  which  necessitates  the  use  of  viscous  circulation 
r(£)  has  indeed  given  better  results,  as  shown  by  Osswald,  K.  Ghia  and  U.  Ghia 
(1992).  Thus,  the  main  emphasis  here  is  on  the  verification  study,  rather  than  on 
validation  of  the  mathematical  model.  In  this  verification  study,  the  accuracy  of 
the  numerical  solutions  obtained  is  examined  critically,  and  the  influence  of  the  dis¬ 
cretized  parameters  of  the  problem  on  the  solution  is  assessed. 

Flow  Configuration  II  described  above  is  used  for  the  verification  study.  The 
parameters  for  this  basic  flow  configuration  II  are  Re  =  45,000,  d+  =  0.2.  Three  dif¬ 
ferent  grid  sizes,  namely  (i)  coarse  grid  with  (330x75)  points,  (ii)  medium  (standard) 
grid  with  (444  x  101)  points,  and  (iii)  fine  grid  (544  x  121)  points  are  used.  The 
minimum  grid  spacings  for  the  standard  grid  in  the  streamwise  and  normal  diretions 
near  the  leading  edge  were  0.54  x  10-3  and  0.26  x  10-4,  respectively.  To  resolve  the 
temporal  scale  of  this  unsteady  separated  flow  problem  and  to  achieve  time  consis¬ 
tency  of  the  solution,  a  value  of  At  =  0.001  was  chosen.  The  code  has  been  fully 
vectoized  and  requires  7.5  micro-seconds  per  time  step,  per  mesh  point,  using  a  single 
processor  of  the  CRAY  Y-MP  8/864  of  The  Ohio  Supercomputer  Center. 

Effect  of  Grid  Stretching  on  the  Solution 

Unique  to  the  present  analysis  is  the  fact  that  the  far-field  boundary  condition  is 
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implemented  at  true  infinity.  As  discussed  earlier  in  the  Section  A. 3.2,  1-D  cluster¬ 
ing  transformations  are  used  to  treat  the  infinite  domain  numerically.  It  is  the  arc 
tangent  transformation  that  permits  the  mapping  of  the  boundary  at  infinity  to  a 
finite  location  in  the  computational  plane.  This  implies  that  the  grids  are  stretched 
exponentially  and  it  is  important  to  examine  the  variation  of  the  dependent  variables 
in  the  far  field  and,  specifically,  their  approach  to  the  imposed  far-field  boundary  con¬ 
ditions,  if  the  solutions  are  to  be  accurate.  Due  to  the  clustering  transformation,  cells 
of  semi-infinite  extent  occur  in  the  physical  plane.  The  effect  of  large  grid  stretch¬ 
ing  on  the  solution  can  be  seen  best  in  the  computational  plane,  and  is  discussed  now. 

Vorticity  attenuates  exponentially  away  from  the  surface  where  it  is  created.  Fig. 
A.6(a-b)  shows  the  vorticity  contours  for  flow  past  a  NACA  0015  airfoil  undergoing 
constant-rate  pitch-up  motion  with  a+  =  0.2,  for  two  different  instantaneous  po¬ 
sitions,  a  =  20.54°  and  36.58°,  respectively.  The  global  views  for  both  instants  are 
shown  in  this  figure,  as  are  the  enlarged  views  in  the  far-field,  which  clearly  show  that 
the  vorticity  approaches  smoothly  to  the  far-field  boundary  condition  of  w/  =  0.  The 
effect  of  grid  stretching  on  the  corresponding  stream  function  is  shown  in  Fig.  A. 7. 
A  close  look  at  Fig.  A.7(a-b)  clearly  shows  that  the  far-field  boundary  is  reached  as 
the  value  of  the  deviational  (i.e.,  the  disturbance)  stream  function  approaches  zero, 
i.e.,  as  rjjf  =  0. 

A  round-off  error  study  was  also  performed,  and  the  procedure  is  as  follows.  In 
the  Dirichlet  Poisson  equation  for  Vs/5,  Eq.  (A. 14),  is  assumed  to  be  known  ana¬ 
lytically  as  any  one  of  the  two  functions  prescribed  in  Fig.  A. 9.  With  this  function, 
the  source  term  involving  w/  is  computed  numerically  using  a  finite-difference  stencil, 
with  extended  precision.  The  computed  values  of  wj  are  used  to  solve  the  stream 
function  equation,  Eq.  (A.  14),  and  the  numerically  computed  values  of  are  com¬ 
pared  with  the  corresponding  known  analytical  values.  The  difference  in  the  two 
solutions  represents  the  round-off  error  on  the  present  (444  x  101)  grid.  Fig.  A. 9 
clearly  shows  that  the  maximum  error  never  exceeds  1 0~  ‘  for  a  solution  field  which 
lies  between  0  and  1.  The  error  also  does  not  exceed  the  maximum  truncation  error 
of  O(10-4).  In  addition,  the  metric  coefficients  used  in  the  calculations  are  shown  in 
Fig.  A. 8;  these  are  well  behaved  as  expected. 

Grid-Refinement  Study 
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The  grid  distribution  for  all  of  the  three  grids  mentioned  earlier,  in  the  proximity 
of  the  NACA  0015  airfoil  are  depicted  in  Fig.  A.  10.  For  the  airfoil  at  zero  angle  of 
attack  in  flow  configuration  II,  the  instantaneous  vorticity  contours  near  the  TE  and, 
subsequently,  in  the  near  wake,  are  presented  here.  They  show  that  the  flow  structure 
consisting  of  nearly  symmetric  vortices  at  the  TE  prevails  for  all  of  the  grids  at  t  = 
3.5.  As  t  increases  to  4.5,  the  results  of  the  finer  grid  show  that  the  vortex  symme¬ 
try  is  broken,  that  is  the  Hopf  bifurcation  has  just  occurred  forming  the  asymmetric 
needle-type  vortical  structures  of  opposite  sign  in  the  flow  field.  At  still  larger  time  t 
=  5.5,  the  flow  field  results  of  medium  grid  show  the  vortex  shedding  process  which  is 
not  only  lagging  the  fine  grid  results  but  it  also  is  out  of  phase  with  the  latter  results. 
The  coarser  results  still  continue  to  show  further  growth  of  the  symmetric  vortical 
structures  near  the  TE.  For  truly  unsteady  flow,  it  is  difficult  to  carry  out  this  grid 
refinement  study,  but  in  light  of  present  results  the  degree  to  which  the  differences 
are  observed  is  clarified. 

Once  the  asymptotic  state  is  achieved  for  the  flow  field  with  the  airfoil  at  zero 
angle  of  attack,  the  constant-rate  pitch-up  motion  with  <i+  =  0.2  was  initiated.  The 
results  in  terms  of  the  instantaneous  vorticity  contours  for  15.96°  <  a  <  35.44° 
are  delineated  in  Fig.  A. 11.  The  coarse-grid  results  are  quantitatively  different.  It 
should  also  be  stated  that,  as  compared  to  the  results  in  Fig.  A.  10,  the  results  for 
the  medium  grid  do  show  departure  from  the  fine-grid  results,  although  qualitatively 
the  two  structures  are  very  similar.  Thus,  the  results  in  Fig.  A.  11  suggest  to  further 
refine  the  grid  and  see  if  the  changes  are  smaller  than  the  previous  refinement  in  which 
grid  size  was  changed  from  (444  x  101)  to  (544  x  121).  Unfortunately,  the  results 
with  the  (544  x  121)  are  computationally  very  costly,  with  the  total  CPU  time  for 
the  fine  grid  case  nearly  5  times  higher  than  for  the  medium  grid  case  using  a  single 
processor  on  the  CRAY  Y-MP  8/864  at  the  Ohio  Supercomputer  Center.  Next,  the 
Cp-distribution  at  the  surface  as  well  as  wall  vorticity  are  plotted  in  Figs.  A.12(a-b) 
for  a  =  14.81°  and  29.71°,  respectively.  For  a  =  14.81°,  where  the  flow  is  mostly 
attached,  the  surface  Cp-distribution  for  all  three  grids  is  in  conformity  and  it  even 
agrees  well  with  the  experimental  data  of  Walker,  Helin  and  Strickland  (1985).  The 
corresponding  wall  vorticity  has  a  similar  behavior,  as  shown  in  Fig.  A. 12(a);  how¬ 
ever,  there  are  no  experimental  data  available  for  wall  vorticity.  On  the  other  hand, 
at  a  =  29.71°,  with  massively  separated  flow  regions,  the  surface  Cp-distributions 
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show  significant  changes,  although  qualitatively  they  do  conform  to  the  experimen¬ 
tal  results  given  by  Walker,  Helin  and  Strickland  (1985).  In  Fig.  A.12(b),  the  wall 
vorticity  on  the  suction  surface  does  show  significant  departures,  again  questioning 
the  grid  independence  of  the  results  between  the  medium  and  fine  grid  at  this  level 
of  the  grid  size.  The  only  satisfactory  comparison  is  that  for  Cl,  as  shown  in  this 
figure.  It  depicts  that  the  coarse-grid  results  compare  well  with  the  experimental 
data  of  Walker,  Helin  and  Strickland  (1985).  However,  the  medium-  and  fine-grid  re¬ 
sults  predict  higher  C^-distribution  in  the  massively  separated  flow  regime.  A  careful 
examination  of  the  experimental  set-up  reveals  that  Walker  et  al.  (1985)  had  only~\ 

36  surface  static  pressure  taps,  as  compared  to  about  207  grid  points  in  the  present  j 
simulation  with  (444  x  101)  grid  points.  Thus,  better  resolution  in  the  computa-  / 
tional  simulation  permits  determination  of  each  of  the  individual  vortices  that  evolve 
and  carries  with  it  somewhat  higher  values  for  (^-distribution  as  compared  to  the 
corresponding  experimental  values.  Hence,  for  theJinestjjndjas^  ^ 

values  of  Cl  are  seen  in  this  figure. 


Comparison  with  Available  Navier-Stokes  Results 


Simulation  results  are  also  obtained  for  a  sinusoidally  pitching  NACA  0012  airfoil 
with  Re  =  5,000,  a  ~  10°(1  —  cos  kt),  and  reduced  frequency  based  on  chord  length,  k 
=  1.0.  Some  experimental  data  is  available  for  this  configuration  from  Werle  (1976). 
Also,  Mehta  (1977)  has  provided  carefully  simulated  Wavier-  Stokes  results  using  the 
NS  equations  in  terms  of  (u/,  0).  He  has  used  an  O-grid  topology  and  provided  detailed 
results.  As  seen  in  Fig.  A.  13,  the  streaklines  for  the  present  results  compare  favorably 
with  the  experimental  data  of  Werle  (1976)  and  show  a  flow  structure  representative 
of  unsteady  flows.  The  flow  is  separated  over  the  entire  suction  surface.  Also  shown 
in  this  figure  are  the  contours  of  instantaneous  vorticity  and  the  velocity  vectors. 
The  better  resolution  of  the  present  results  is  clear  and  these  results  conform  well 
with  those  of  Mehta  (1977).  For  this  sinusoidally  pitching  airfoil,  the  instantaneous 
stream-function  contours  are  compared  in  Fig.  A.  14.  The  results  first  show  one 
instant  for  a  =  18.59°  in  the  upward  stroke,  which  goes  up  to  a  =  20°,  and  then  four 
values  of  a  for  the  return  stroke.  For  the  stream  function,  which  is  not  very  sensitive, 
the  agreement  is  good.  For  the  same  values  of  a,  comparison  is  also  provided  for  the 
instantaneous  vorticity  contours  in  Fig.  A.  15.  The  qualitative  agreement  between 
the  two  sets  of  results  is  very  good,  although  the  better  resolution  of  the  present 
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results  due  to  the  use  of  a  finer  grid,  made  possible  by  a  better  and  larger  computer 
is  very  evident.  The  complete  evolution  of  the  dynamic  stall  event  is  vivid,  and  some 
of  the  interactions  that  take  place  can  be  inferred  from  this  figure.  Finally,  the  C i- 
distribution,  as  well  as  the  surface  pressure  coefficients  for  two  values  of  a,  namely, 
a  =  18.57°  during  pitch-up  and  a  =  11.11°  during  pitch  down,  are  shown  in  Fig.  A.  16. 
Also  shown  in  this  figure  are  the  results  of  Sankar  and  Tassa  (1980).  In  general,  Cl 
compares  qualitatively  with  the  existing  results,  but  there  are  significant  deviations. 
For  the  time  being,  this  completes  the  verification  study;  additional  verification  with 
experimental  data  is  provided  in  the  next  sub-section. 

Influence  of  Initial  State  on  the  Flow  Field 

For  flow  configuration  II,  results  are  initially  obtained  to  determine  the  asymptotic 
state  at  fixed  angle  of  attack.  Since  this  asymptotic  state  is  not  unique,  the  effect  of 
the  initial  state  of  flow  will  prevail  in  the  final  solution.  To  determine  this  effect,  three 
different  asymptotic  states  are  selected,  and  correspond  to  states  (b),  (c)  and  (d)  on 
the  Cx-history  curve  in  Fig.  A.17  corresponding  to  a  =  0°.  With  these  states  aa 
the  starting  solutions,  three  different  calculations  are  made  for  constant-rate  pitch-up 
motion  and  the  results  are  compared  in  Fig.  A.17.  The  lift  coefficient  distribution,  as 
well  as  the  viscous  circulation,  for  all  three  cases  do  conform  with  each  other  and  do 
not  show  strong  influence  of  the  initial  state.  The  instantaneous  vortidty  contours  are 
also  shown  in  this  figure  and,  although  there  is  a  phase  difference  in  the  development 
of  the  wake,  the  overall  influence  is  minimal. 

Flow  Structure  for  Configuration  EE:  Re  —  45,000,  No  Suction 

For  this  configuration,  Walker,  Helin  and  Strickland  (1985)  have  provided  flow  visu¬ 
alization  photographs  for  six  selected  values  for  angle  of  attack.  These  photographs 
very  vividly  show  initially  the  region  of  separated  flow  near  LE  and  TE  and,  sub¬ 
sequently,  the  formation  of  a  highly  energetic  dynamic  stall  vortex.  Streaklines  are 
plotted  in  Fig.  A.  18  at  a  =  24°, 27°, 30°, 40°  and  47°  and  they  compare  well  with  the 
flow  visualization  data.  Further,  the  C^-history,  as  well  as  the  surface  distribution 
of  Cp  are  compared  with  the  experimental  data  in  Fig.  A.19(a,b).  As  seen  in  Fig. 
A.  19(a),  the  Cx-distribution  agrees  qualitatively,  but  shows  higher  values  of  Cl-  It 
was  also  pointed  out  that  better  agreement  with  experimental  data  is  obtained  for  the 
coarse-grid  simulation  using  (330  x  75)  points.  Also  shown  here  is  Cx-distribution 
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predicted  by  Visbal  and  Shang  (1989);  this  too  agrees  better  with  the  experimen¬ 
tal  data.  In  the  opinion  of  the  present  authors,  this  is  again  due  to  their  relatively 
somewhat  coarser  grid  of  (203  x  101)  points,  coupled  with  their  treatment  of  the 
far-field  boundary  condition.  Next,  the  coefficient  of  surface  pressure  distribution  is 
compared  for  two  values,  a  =  14.81°  and  29.71°,  in  Fig.  A.19(b).  As  seen  here, 
the  comparison  with  the  experimental  data  of  Walker,  Helin  and  Strickland  (1985) 
is  good  for  a  =  14.81°,  whereas,  for  a  =  29.71°,  the  two-distributions  depart  signifi¬ 
cantly  at  the  suction  surface.  It  is  again  the  opinion  of  the  authors  that  this  is  due  to 
three-dimensionality  in  the  flow  field.  The  results  presented  so  far,  provide  additional 
verification  of  the  present  analysis. 

Next,  the  detailed  flow  structure  is  depicted  in  Figs.  A.20-A.21  using  the  instan¬ 
taneous  vorticity  contours  with  a  varying  from  11.37°  to  29.71°.  Here,  global  flow 
structure  information  is  made  available,  along  with  an  enlarged  view  showing  the  LE 
flow  structure.  As  the  airfoil  is  pitched-up,  flow  reversal  takes  place  near  the  LE  as 
shown  in  this  figure.  Subsequently,  the  flow  reversal  evolves  into  secondary,  tertiary 
and  quaternary  vortices.  The  eruption  of  this  secondary  vortex  structure  initiates 
the  formation  of  the  energetic  dynamic-stall  vortex.  K.  Ghia,  Yang,  Osswald  and  U. 
Ghia  (1992)  have  provided  further  details  of  this  behavior. 

Flow  Structure  for  Configuration  III:  Re  =  45,000  with  Modulated  Suc¬ 
tion/Injection 

In  an  attempt  to  control  the  formation  of  the  dynamic  stall  vortex,  two  specific 
MSI  controls  were  used.  These  are  Case  1  with  Ui/f/oo  =  0.045  and  Case  2  with 
Vj/f/oo  =  0.035,  as  presented  in  Fig.  A. 22.  The  ramp  functions  used  for  these  two 
configurations  to  bring  the  magnitude  of  suction  velocity  to  its  final  magnitude  are 
shown  in  this  figure;  these  magnitudes  are  held  constant  thereafter.  Also  shown  in 
this  figure  is  the  viscous  circulation  for  configuration  III  as  well  as  for  both  cases  of 
configurations  II.  The  start  of  the  MSI  is  also  indicated.  In  addition,  the  velocity 
vectors  are  shown  for  both  cases  for  three  different  values,  a  =  23.98°,  25.12°  and 
26.27°.  Also  shown  here  are  the  velocity  vectors  as  well  as  Cp-distribution  on  the 
suction  surface  for  flow  without  and  with  MSI.  In  Fig.  A. 23,  the  instantaneous 
vorticity  contours  as  well  as  the  wall  vorticity  are  shown  for  configurations  II  and  III. 
In  this  figure,  results  are  also  plotted  for  one  additional  value  of  a  =  22.83°.  The 
instantaneous  vorticity  contours  for  various  a  show  that  the  dynamic  stall  vortex  does 
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not  develop  with  MSI  activated.  Further,  the  controlled  wall  vorticity  does  not  show 
sharp  spikes  leading  to  suppresion  of  the  formation  of  the  stall  vortex.  The  results 
showing  the  effect  of  MSI  on  Cl,  Cq,  Cm  and  L/D  are  most  significant,  since  they  show 
directly  whether  or  not  the  desired  control  is  achieved.  Fig.  A. 24(a)  shows  that  there 
is  not  much  significant  change  in  Cl  ,  but  as  seen  in  Fig.  A. 24(b),  Co  is  considerably 
reduced  by  MSI  for  Case  1  at  higher  a.  Although  Cm  distribution  does  not  show  a 
clear  trend,  in  Fig.  A. 24(c),  L/D  in  Fig.  A. 24(d)  increases  rapidly  as  compared  to 
that  for  the  flow  without  MSI.  Finally,  controlled  modulated  suction/injection  is  also 
attempted  using  configuration  III,  with  a  more  gradual  application  of  MSI  control 
as  shown  in  Fig.  A. 25.  Compared  to  Case  1,  the  viscous  circulation  is  sufficiently 
reduced,  for  this  new  Case  3.  As  before,  the  velocity  vectors  as  well  as  Cp-distribution 
on  the  suction  surface  are  also  plotted.  The  corresponding  instantaneous  vorticity 
contours  and  wall  vorticity  are  shown  in  Fig.  A. 26.  The  results  for  Case  3  with  its 
gradually  implemented  MSI  are  very  similar  to  those  for  Case  1.  Finally,  the  effect 
of  gradually  applying  MSI  does  not  show  any  significant  changes  as  compared  to  the 
standard  Case  1. 

2.A.6  Conclusion 

An  unsteady  NS  analysis  is  developed  using  the  (u>/,  V?  ,  T(t))  formulation.  Inclusion 
of  viscous  circulation  permits  accurate  treatment  of  the  far-field  boundary  condi¬ 
tion.  The  implicit  ADI-BGE  numerical  technique  provides  a  uniformly  second-order 
accurate  solution,  with  higher  accuracy  for  the  nonlinear  convective  terms.  The  grid- 
independence  study  showed  that,  for  this  constant-rate  pitch-up  motion  in  the  regimes 
where  the  flow  is  massively  separated,  the  finest-grid  solutions  still  show  more  than 
desired  deviation  from  those  of  the  medium  grid  used  here  to  obtain  all  the  solutions. 
Also,  since  the  experimental  pressure  distribution  is  measured  using  a  limited  number 
of  pressure  taps,  the  coarse-grid  CL-distribution  agrees  well  with  these  results.  The 
role  of  unsteady  separation  was  discussed  by  the  authors  in  detail  in  Ref.  [A.  12]  and 
is  therefore  not  described  here.  Based  on  that  earlier  description,  a  modulated  suc¬ 
tion/injection  control  strategy  was  developed  and  used  in  this  study.  The  suppression 
of  LE-flow  separation  provides  a  very  effective  means  of  controlling  the  formation  of 
the  dynamic  stall  vortex,  even  with  small  volumetric  suction  rate  S  =  -0.00237. 

The  authors  believe  that  the  2-D  analysis  developed  here  is  accurate  and  that  any 
further  discrepancies  observed  are  due  to  three-dimensionality  effect.  An  effort  should 


26 


be  made  to  accurately  simulate  the  three-  dimensional  flow  past  this  2-D  geometry. 
Effort  is  also  necessary  to  incorporate  an  adaptive-grid  refinement  in  this  analysis,  in 
order  to  obtain  higher-Re  solutions  and  also  resolve  the  time  scale  of  the  eruption  of 
the  secondary  structure  of  0 (Re~2^11).  The  three-dimensional  analysis  should  also  be 
coupled  with  a  large-eddy  simulation  technique  in  order  to  obtain  higher-Re  solutions 
of  greater  practical  interest. 
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2.B  Three-Dimensional  Iterative  Technique  with 
Multi-Grid  Acceleration 

2.B.1  Introduction 

Three-dimensional  unsteady  incompressible  flows  are  of  general  practical  and  compu¬ 
tational  interest.  Modern  computers  have  enabled  simulation  of  more  complex  flows 
occurring  in  engineering  applications.  A  number  of  numerical  algorithms  have  been 
developed  to  solve  incompressible  3-D  Navier-Stokes  equations,  using  either  primitive 
variables  or  vorticity-based  formulations  (Refs.  [B.6,  11,  21,  22,  23,  24]). 

Efforts  are  regularly  pursued  for  the  optimization  of  numerical  simulation  tech¬ 
niques  to  obtain  enhanced  computational  accuracy  and  efficiency.  The  present  study, 
thus,  is  aimed  at  developing  an  efficient  algorithm  for  the  velocity  problem  emerg¬ 
ing  from  the  vorticity-velocity  formulation  of  3-D  Navier-Stokes  equations,  using  fine 
grids.  Multigrid  (MG)  techniques  have  provided  an  efficient  methodology  to  achieve 
this.  The  present  sub-section  describes  the  development  of  an  MG  method  and  its 
application  in  unsteady  three-dimensional  incompressible  flow  simulation  using  the 
vorticity-velocity  (Q-V)  formulation  of  the  Navier-Stokes  equations.  Second-order 
accuracy,  both  in  time  and  space,  is  obtained  via  central-differencing  the  vortidty- 
transport  equations  and  the  elliptic  velocity  problem. 

2.B.2  Vorticity- Velocity  Formulation  of  3-D  Navier-Stokes 

Equations 

The  vorticity-velocity  (w-V)  formulation,  as  opposed  to  the  primitive-variable  [V -p) 
formulation,  is  selected  in  the  present  study  because  the  former  has  the  impressive  fea¬ 
ture  that  the  spin  dynamics  of  a  fluid  particle  (represented  by  the  vorticity-transport 
equations)  is  naturally  separated  from  the  kinematic  motion  of  the  fluid  particle  (rep¬ 
resented  by  the  elliptic  velocity  problem).  It  is  recognized  that  the  vorticity  is  a  linear 
function  of  velocity,  whereas  the  pressure  is  a  nonlinear  function  of  velocity.  There¬ 
fore,  the  resultant  numerical  methods  are  easier  to  deal  with,  from  a  computational 
point  of  view.  Furthermore,  with  appropriate  interpretation  of  the  variables,  the 
vorticity-velocity  formulation  retains  its  form  even  in  non-inertial  coordinate  systems. 

As  discussed  in  Refs.  [B.13,  18,  19],  the  non-dimensionalized  incompressible 
Navier-Stokes  equations  can  be  written  as  follows  in  terms  of  vorticity-velocity  vari- 
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(B.l) 


ables  in  general  vector  conservation  form: 

du  -  1 

—  +  V  x  (w  x  V)  =  *V  x  (V  x  u>) 
at  rie 

where  R z—LrVrIv  is  the  Reynolds  number,  based  on  an  appropriate  reference  length 
Lr  and  velocity  Vr. 

The  kinematic  velocity  problem  (KVP)  is  formulated  directly  from  the  continuity 
equation  and  the  definition  of  the  vorticity  vector: 

V*?  =  0,  (B.2) 

Vxi^  =  u  (B.3; 

subject  to  constraints: 

JjV'ndS  =  0  (B.4) 

and 

V.JsO.  (B.5) 

Eqs.  (B.4)  and  (B.5)  are  necessary  and  sufficient  conditions  for  the  existence  of  a 
unique  solution  for  Eqs.  (B.2)  and  (B.3).  Consequently,  the  3-D  velocity  problem 
is  completely  described  and  determined  by  Eqs.  (B.2-B.5)  and  the  boundary  condi¬ 
tions.  It  is  important  to  note  that  the  constraint  Eq.  (B.5)  is  time- invariant,  whereas 
applying  the  divergence  operator  on  Eq.  (B.l)  yields  that  <9(V  «u ))jdt  =  0;  together, 
these  ensure  that  the  velocity  is  determined  uniquely  for  t  >  0. 

The  governing  equations  and  constraints,  Eqs.  (B.1-B.5),  are  then  expressed  in  a 
general  curvilinear  orthogonal  coordinate  system  (C^2A3)-  Covariant  components 
of  both  velocity  and  vorticity  are  employed.  This  permits  the  strong  conservation-law 
form  of  the  governing  equations  to  be  maintained  and  the  constraints  to  be  satisfied 
analytically  in  the  generalized  coordinate  system.  Another  significant  consequence 
of  the  covariant  decomposition  of  both  velocity  and  vorticity  vectors  is  that  the  gen¬ 
eralized  form  of  the  curl  Eq.  (B.3)  contains  no  variable  coefficients.  Furthermore, 
the  components  of  the  vorticity-transport  equation  are  characterized  by  the  fact  that 
mixed  derivatives  appearing  in  the  equation  for  one  component  (e.g.,  a-vequation) 
involve  only  the  other  two  components  of  vorticity  (u2  and  u>3);  this  allows  implicit 
treatment  of  these  terms  within  an  approximately-factored  solution  scheme  (Refs. 
[B.13,  18,  19]). 


29 


2.B.3  Numerical  Scheme 


Second-order  finite  differencing  is  used  in  the  discretization  of  both  vorticity-transport 
equations  and  the  kinematic  velocity  problem  on  a  staggered  grid.  The  velocity  com¬ 
ponents  are  defined  at  the  centroid  of,  and  perpendicular  to,  the  cell  surfaces.  The 
vortidty  components  are  tangent  to  the  grid  lines  and  are  located  at  the  mid-point  of 
the  cell  edges.  This  particular  staggered  arrangement  has  the  advantage  that  a  set  of 
physically  consistent  algebraic  equations  may  be  derived,  that  satisfy  the  constraint 
Eqs.  (B.4)  and  (B.5). 

First  of  all,  The  boundary  conditions  on  velocity  must  be  such  that  the  global 
mass  is  conserved  (Eq.  (B.4j).  Strong  conservation-law  form  of  Eq.  (B.2)  guarantees 
that  the  integral  constraint  (Eq.  (B.4))  will  be  satisfied  if  the  equation  of  continuity 
at  each  computational  cell  is  satisfied. 

Secondly,  the  discrete  counterpart  of  Eq.  (B.5)  will  be  satisfied  at  any  time  during 
the  evolution  of  the  flow  if  it  is  satisfied  initially  and  the  Algebraic  vorticity-transport 
equation  is  solved  by  a  direct  method.  However,  in  the  present  work,  the  vorticity- 
transport  equations  are  solve-  /  a  modified  ADI  method  developed  by  Osswald  et 
al.  (1988).  The  approximate  factorization  creates  a  splitting  error  of  order  0(At3); 
the  time-marching  solution  of  the  vorticity  field  will  then  be  divergence-free  up  to  the 
same  order.  Consequently,  the  convergence  of  the  iterative  solution  for  velocity  field 
will  be  limited  by  this  error. 

For  the  velocity  problem,  Osswald  et  al.  (1987)  had  employed  a  direct  solver,  Block 
Gaussian  Elimination  (BGE),  to  compute  the  velocity  field.  Due  to  the  limitation  of 
the  current  core  memory  of  the  available  supercomputers,  application  of  the  direct 
solver  is  restricted  by  the  allowable  grid  size.  Therefore,  iterative  schemes,  widely 
utilized  in  the  solutions  of  PDE’s  (Refs.  [B.l,  3,  9,  10,  17])  are  explored  in  the  present 
study.  Table  B.l  compares  the  requirements  of  storage  between  the  direct  solver 
BGE  and  the  Multi-Grid  Distributive  Gauss-Seidel  (MG-DGS)  iterative  algorithm 
developed  in  the  present  study  for  the  velocity  problem,  it  is  evident  that  the  iterative 
method  will  ease  the  memory  demands.  However,  iterative  methods  frequently  exhibit 
slow  convergence  rates. 
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Multigrid  Method 

Multigrid  methods  can  help  to  accelerate  the  convergence  of  the  iterative  process 
for  the  discretized  system  of  Eqs.  (B.2-B.3),  because  the  system  is  of  elliptic  type 
Refs.([B.3,  B.8]).  it  is  known  that  any  iterative  scheme  is  only  effective  in  its  ini¬ 
tial  stage,  when  oscillatory  (high  frequency)  error  components  are  dominant.  The 
slowdown  thereafter  is  caused  by  the  ineffectiveness  of  the  relaxation  to  smooth  (low 
frequency)  errors  on  the  current  fine-grid.  By  switching  to  a  coarser  grid  level,  these 
low-frequency  error  components  will  become  oscillatory  on  the  coarse  grid,  and  will 
be  smoothed  out  efficiently  by  the  relaxation. 

Although  variable  coefficients  are  possible  due  to  coordinate  transformation  in 

Eq.  (B.2),  the  system  of  Eqs.  (B.2-B.3)  is  always  linear.  As  a  result,  the  Correction 

Scheme  (CS)  of  Brandt  (1979,  1984)  is  employed  in  the  current  study.  A  stable  and 

efficient  Distributive  Gauss- Seidel  relaxation  scheme  is  developed  to  smooth  out  os- 

£ 

cillatory  errors.  Only  the  residuals  on  the  find  grid  need  be  restricted  to  the  coarse 
grid.  The  correction  is  obtained  by  relaxation  on  the  coarse  grid;  it  is  then  interpo¬ 
lated  (prolongated)  back  to  the  fine  grid  to  correct  the  old  fine  grid  approximation. 
The  restriction,  interpolation  and  coarse-grid  operators  will  be  discussed  later  in  the 
sub-section. 

Standard  full  coarsening  is  used  to  define  coarser  grids.  As  a  result,  one  coarse 
computational  cell  contains  eight  fine  cells.  Because  of  the  staggered  grid  arrange¬ 
ment,  the  coarse-grid  function  locations  never  coincide  with  the  locations  of  the  fine- 
grid  functions.  The  Multigrid  cycles  used  in  the  present  study  are  interacted  adap¬ 
tively  (Accommodative  Cycle).  Since  the  flow  in  the  present  study  is  time  dependent, 
the  iterative  solution  of  the  velocity  field  is  repeated  at  each  temporal  step.  The  so¬ 
lution  at  the  previous  time  step  serves  as  a  good  approximation  for  the  solution  at 
the  new  time  step.  Thus,  the  MG  cycle  starts  at  the  finest  grid,  visiting  the  coarser 
grids  successively  to  obtain  the  correction. 

Relaxation  Scheme 

To  arrive  at  a  fast  MG  algorithm,  the  discretization  of  the  boundary- value  problem 
should  be  stable  in  the  first  place.  The  central  differencing  scheme  on  the  staggered 
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grid  used  in  the  present  study  for  the  velocity  leads  to  finite-differenced  equations 
that  are  stable  and  elliptic  (Ref.  [B.3j). 


The  standard  pomtwise  Gauss-Seidel  (GS)  relaxation  is  usually  efficient  for  a 
single  equation.  However,  the  discretized  system  for  the  velocity  problem  does  not 
contain  a  one-to-one  correspondence  between  the  unknowns  and  the  equations.  The 
equations  are  coupled  to  each  other,  and  therefore,  a  pointwise  GS  scheme  would  not 
work  effectively.  An  alternate  scheme,  Distributive  Gauss-Seidel  (DGS)  scheme,  is 
thus  recommended  for  systems  of  equations  (Ref.  [B.l]).  In  the  DGS  scheme,  each 
discretized  equation  in  the  system  is  satisfied  in  turn  by  distributing  its  residual  to 
all  the  unknown  variables  associated  with  this  equation.  This  is  done  in  such  a  way 
that  the  residuals  of  the  remaining  equations  are  kept  unaltered.  Mathematically  , 
the  system  for  Eqs.  (B.2-B.3),  can  be  written  symbolically  as  a  matrix  operator  form 
LV  =  /  with  L  =  (V#,  Vx)r.  It  is  then  reformulated  as  follows: 

LLw  =  /,  (B.6) 

where  the  operator  LL,  preferably,  has  a  diagonal  (or  block  diagonal)  structure.  The 
Gauss-Seidel  iteration  is  then  formulated  with  respect  to  the  new  unknown  w  (called 
“ghost  unknown”  by  Brandt  (1984).  A  new  iteration  V  —  Lw  will  yield  the  original- 
grid  function  V.  In  fact,  the  new  unknown  w  is  neither  stored  nor  calculated  in  the 
actual  program. 


The  advantage  that  the  DGS  scheme  gains  is  that  LL  usually  possesses  better 
smoothing  rate  than  L  itself.  For  the  present  velocity  problem,  the  operator  LL  will 
become 


LL  = 


V2  0 
0  V  x  Vx 


(B.7) 


if  L  is  chosen  as  (V,  Vx).  It  is  seen  that  the  operator  LL  contains  the  second-order 
Laplacian;  better  smoothing  rate  is  expected  when  relaxing  Eq.  (B.6). 


Specifically,  the  continuity  equation  is  scanned  and  satisfied,  while  keeping  un¬ 
changed  the  residuals  of  all  neighboring  curl  equations.  The  distributed  amount 
should  be  suitably  modified  for  boundary  cells  where  only  interior  velocities  are  up¬ 
dated  to  satisfy  the  continuity  equation  for  the  boundary  cells.  Red-Black  ordering 
is  employed  for  the  relaxation  of  the  continuity  equation  to  achieve  vectorization. 
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For  the  curl  equations,  the  present  DGS  scheme  groups  all  three  curl  equations 
and  relaxes  them  simultaneously  in  the  prescribed  order,  keeping  all  the  residuals  of 
the  neighboring  continuity  equations  unchanged.  Although  the  scheme  turns  out  to 
be  a  little  cumbersome  algebraically,  the  principle  is  exactly  the  same  as  for  relaxing 
the  continuity  equation.  Four-color  ordering  is  used  to  eliminate  data  dependence, 
and  achieve  vectorization  of  the  relaxation  sweep  for  all  three  curl  equations. 

Lexicographic  ordering  is  also  used  in  the  present  study  as  well  as  red-black/four- 
color  ordering.  It  is  found  that  the  smoothing  rates  for  both  orderings  differ  very 
little  for  the  present  velocity  problem. 

Interpolation 


The  operator  interpolates  the  coarse-grid  ( M  —  1)  function  to  the  fine  grid 
(M).  The  order  of  the  interpolation  should  be  dictated  by  the  order  of  the  difference 
(or  differential)  equations  (Refs.  [B.l,  B.8]).  In  addition,  the  integral  constraint, 
Eq.  (B.4),  must  be  satisfied  on  all  grid  levels  so  that  consistency  is  maintained. 
Trilinear  interpolation  for  the  interior  grid  points  and  linear  interpolation  for  the 
near-boundary  points  are  constructed  to  fulfill  these  requirements.  Each  of  the  three 
velocity  components  is  interpolated  separately.  Then,  the  interpolation  for  VI  follows: 

+  3F1"-A  +  3Kl"r+‘,  +  n"-A+1),  (interior) 

(boundary) 

IS- ,n«rl  =  ^"r‘  (corner)  (B.8) 


for  the  fine-grid  variables  on  the  coarse-grid  plane,  and 

jM  i/iAf-l  _  ^  (  tM  t/ i  A/ —  1  i  rrVf  t/iA/-1  \ 

lM-lv  Li+ljk  -  L,jk  +  ‘Af —1  ’  H+2 jk) 


(B.9) 


for  the  fine-grid  points  on  the  plane  between  two  coarse-grid  planes.  Similar  expres¬ 
sions  are  developed  for  the  V2  and  V3  components.  In  case  of  variable  coefficients 
as  a  result  of  non-uniform  grid,  the  variable  V\  in  Eqs.  (B.8-B.9)  should  be  taken  as 
the  product  of  the  Vx  and  the  corresponding  coefficient.  Such  a  scheme  ensures  that 
the  integral  constraint  is  satisfied  on  all  grid  levels. 


Restriction  jT^"1 
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In  the  Correction  Scheme,  only  residuals  need  to  be  transferred  to  coarser  grids. 
Full  weighting,  which  averages  all  of  the  neighboring  fine-grid  residuals,  is  especially 
important  for  difference  equations  with  varying  coefficients.  This  is  exactly  the  case 
for  the  velocity  problem  in  the  present  study.  Thus,  full  weighting  (FW)  residual 
transfer  is  constructed  for  all  the  residuals  in  the  study. 

As  the  finite-difference  equations  for  the  velocity  problem  are  centered  at  different 
locations,  residuals  are  transfered  separately  for  each  equation.  Specifically,  for  the 
continuity  equation,  the  restriction  operator  is: 

-  xfVJ'xJ'MSl,  (B.10) 

where  =  j(/"  ijk  +  /i+ij  J  is  t^e  mid-point  average  operator  for  any  grid  func¬ 

tion  f-fk  at  level  M,  and  rd  is  the  residual  of  the  continuity  equation.  The  coarse- grid 
source  term  for  the  continuity  equation  is  thus  the  average  of  the  neighboring  eight 
fine-grid  residuals. 

Similarly,  the  full  weighting  for  the  residual  transfer  of  W\  (discretized  o?i)  com¬ 
ponent  of  the  curl  equation  is  written  symbolically  as: 

\rW\)ijk  —  Pj  Pk  Pj  Mfc  Mi  VWIJijk-  (B-li) 

This  says  that  the  coarse-grid  source  term  W iff*.-1  is  the  weighted  average  of  the 
eighteen  neighboring  fine-grid  residuals.  The  residuals  for  W2  and  IT 3  components 
of  the  curl  equation  are  treated  in  a  similar  manner. 

The  above  restriction  operator  also  preserves  the  two  compatibility  conditions, 
given  by  Eqs.  (B.4)  and  (B.5),  throughout  all  the  coarse-grid  levels,  as  long  as  they 
are  met  on  the  finest  grid. 

Coarse-Grid  Operator  LM~l 

The  guideline  in  constructing  the  coarse-grid  matrix  operator  LM~  1  is  that  LM'1 
should  be  a  proper  homogenization  of  the  fine-grid  operator  LM .  One  way  to  achieve 
this  is  to  define  LM~l  in  the  same  manner  as  LM  from  the  corresponding  differential 
operators.  This  approach  is  employed  in  the  present  study.  Therefore,  no  additional 
computations  are  needed  for  defining  Lw_1,  except  when  varying  coefficients  occur 


34 


in  Eq.  (B.2)  because  of  coordinate  transformation.  In  this  case,  some  appropriate 
averaging  is  needed  to  define  the  coarse-grid  coefficients.  Since  the  varying  coefficients 
for  the  fine-grid  are  due  to  the  coordinate  transformation,  the  weighted  average  of  the 
fine-grid  coefficients  is  naturally  a  proper  definition  for  the  coarse-grid  coefficients, 
and  thus  for  the  coarse-grid  operator  LM~X.  The  so-called  minimal  weighting  (Ref. 
[B.3])  is  employed  for  forming  the  coarse-grid  coefficients,  as  for  example, 

cu&-‘  -  Mj'dJ'cnj?,  (B.i2) 

where  Gllijk  is  the  discretized  metric  coefficient  occurring  in  Eq.  (B.2)  due  to  coor¬ 
dinate  transformation. 

2.B.4  Model  Velocity  Problem 

Numerical  experiments  have  been  conducted  using  the  DGS  relaxation  scheme  de¬ 
veloped  together  with  the  MG  procedure.  A  test  function  cos(sin(xyz ))  is  used  to 
generate  a  ‘velocity’  field  in  the  Cartesian  domain.  The  divergence  and  curl  of  this 
vector  field  are  then  computed  using  the  discretized  form  of  the  governing  equations 
(Eqs.  (B.2-B.3))  and  provide  the  appropriate  source  terms  for  these  equations.  Then, 
the  multigrid  method  developed  is  used  to  solve  these  equations  with  modified  source 
terms  and  determine  ‘velocity’  field,  starting  with  an  appropriate  initialization  (zero 
is  taken  in  the  present  test  case),  until  the  prescribed  convergence  criterion  is  satis¬ 
fied.  Both  uniform  and  non-uniform  grids  have  been  used.  The  non-uniform  grids 
were  generated  by  algebraic  clustering. 

Figures  B.1-B.3  demonstrate  graphically  how  the  MG-DGS  scheme  improves  the 
convergence  rate,  relative  to  the  single-grid  (SG)  results,  for  both  uniform  and  clus¬ 
tered  (25  x  25  x  25)  grids.  Table  B.2  gives  some  numerical  values  for  the  MG-DGS 
scheme  for  several  uniform  and  non-uniform  grids. 

It  is  seen  that  the  CPU  time  increases  linearly  and  the  convergence  rate  is  main¬ 
tained  at  about  QA2/WU  for  uniform  grids  with  the  size  increasing  from  (25  x  25  x  25) 
to  (33  x  33  x  33)  to  (49  x  49  x  49).  Similar  results  are  obtained  for  non-uniform 
(25  x  25  x  25)  and  (33  x  33  x  33)  grids  with  moderate  clustering.  The  expected 
multigrid  convergence  rate  is  obtained.  The  convergence  rate  jj,  in  Table  B.2  is  the 
average  of  the  total  rate  defined  by: 

_  /  II  r  II  final  U fWU 
^  '||  f  !| initial 
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where  the  norm  |j  •  ||  is  defined  as  the  root-mean-square  of  the  residuals  and  WU  is 
the  number  of  work  units  executed. 

It  is  seen  that  when  varying  coefficients  occur,  as  a  result  of  grid  clustering,  the 
asymptotic  convergence  rate  degrades  slightly.  Grid  spacings  for  non-uniform  grids 
used  in  Table  B.2  have  differences  of  one  order  of  magnitude.  The  deterioration  is 
greater  for  the  case  with  more  highly  varying  coefficients.  This  is  demonstrated  in 
Fig.  B.3,  where  the  minimal  spacing  is  of  order  (9(1Q~3)  and  maximal  spacing  is  of 
O(K)'1);  i.e.,  a  difference  of  two  orders  of  magnitude  exists  in  the  coefficients.  Further 
numerical  experiments  with  different  test  functions  indicate  that  the  convergence  also 
depends  upon  the  nature  of  the  velocity  field  itself.  It  is  noted  that  the  clustering 
is  strong  near  the  boundary  (wall)  of  the  domain  because  high  gradients  of  the  flow 
variables  (vorticity  and  velocity)  are  expected  here  in  the  actual  viscous  flows.  Figure 
B.4  shows  the  convergence  history  obtained  from  the  actual  flow  simulation  for  the 
3-D  shear-driven  cavity  of  aspect  ratio  unity  using  the  same  clustering  as  in  Fig.  B.3. 
Desirable  convergence  rate  is  observed  for  this  strongly  clustered  grid  (two  orders  of 
magnitude  difference),  thereby  demonstrating  the  robustness  of  the  MG-DGS  scheme. 

2.B.5  Some  Physical  Aspects  of  Unsteady  3-D  Shear. Driven 
Flow 


Shear-driven  flow  in  a  cavity  has  drawn  much  attention  in  the  fluid  dynamics  commu¬ 
nity  because  of  its  richness  of  flow  structure  and  complexity  of  dynamical  phenomena, 
as  well  as  the  simplicity  of  its  geometry  and  well-defined  boundary  conditions.  Care¬ 
ful  experiments  have  been  conducted  by  Koseff  and  Street  (1984),  among  others,  to 
study  the  characteristics  of  this  flow.  The  cavity  has  also  been  chosen  as  a  benchmark 
test  by  many  CFD  practitioners:  for  example,  Refs.  [B.7]  and  [B.21],  for  both  2-D 
and  3-D  flow  simulations. 


The  shear-driven  cavity  flow  has  many  complex  flow  structures.  Multiple  recir- 
culatioi 


^xist  as  a  result  of  strong  viscous  effects  exerted  by  the  walls  and  the 
sharp  corners  in  the  geometry.  Of  special  interest  are  the  three  dimensionality  and 
th£ longitudinal  vortices.  Efforts  are  made  in  the  present  study  to  investigate  theJ3-D 
cavity  flow  using  the  MG-DGS  method  to  validate  the  method  as  well  as  to  study 
the  flow  itself. 
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Figure  B.5  depicts  the  configuration  of  the  3-D  shear-driven  cavity  flow.  The  up¬ 
per  wall  ABFE  is  impulsively  set  in  motion  in  the  direction  parallel  to  itself.  The 
flow  is  thus  symmetric  with  respect  to  the  mid-span  plane.  The  typical  time-averaged 
flow  pattern  in  the  symmetry  plane  is  shown,  schematically,  in  Fig.  B.6.  The  flow 
exhibits  three  significant  recirculation  regions,  in  addition  to  the  primary  vortex:  the 
downstream  secondary  vortex,  the  upstream  secondary  vortex  and  the  upper-surface 
secondary  vortex.  These  features  make  this  flow  a  benchmark  model  for  testing  nu¬ 
merical  methods  for  incompressible  viscous  flows. 

The  formation  of  these  recirculation  regions  is  the  result  of  the  wall  friction  and 
the  stagnation  induced  by  the  corner.  Moreover,  the  separation  zone  downstream 
provides  a  concave  boundary  for  the  primary  flow  coming  down  from  above.  This 
suggests  the  occurrence  of  the  Taylor- Gortler-like  (TGL)  vortex  above  the  separation 
surface.  The  phenomenon  was,  indeed,  observed  in  the  experimental  study  by  Koseff 
et  al.  (1983,  1984).  These  vortices  are  susceptible  to  instability;  they  may  eventually 
lead  to  transition  at  higher  Reynolds  number. 

Examination  of  TGL  Vortices 

In  an  attempt  to  capture  the  TGL  vortices,  a  case  with  higher  Reynolds  number 
Re  =  3300  is  investigated.  The  span  wise  aspect  ratio  is  3:1.  This  corresponds  to  the 
experimental  setup  used  in  Ref.  [B.14],  Figure  B.7  shows  the  visualization  results 
obtained  by  Koseff  et  al.  (1983).  The  TGL  vortices  in  the  lateral  plane  at  x1  =  0.8 
are  clearly  seen  in  Fig.  B.7(b). 

To  resolve  the  viscous  boundary  layer  near  the  wall  and  the  shear  layer  at  the 
upper  surface,  a  non-uniform  (65  x  65  x  33)  grid  is  generated,  with  clustering  in  the 
x1  and  x2  directions.  The  total  core  storage  requirement  for  this  grid  size,  using  the 
ADI+MG-DGS  method,  is  about  5  megawords.  A  time  step  of  A t  =  0.01  is  used  in 
the  marching  solution  for  the  unsteady  flow.  The  convergence  criterion  for  the  finest 
grid  for  the  MG-DGS  method  for  the  velocity  problem  is  taken  to  be  5  x  10-7.  Table 
3  lists  the  CPU  index  for  the  flow  simulation;  comparison  is  also  made  with  the  BGE 
method  (Ref.  [B.19])  for  the  (25  x  25  x  25)  uniform  grid  at  Re  =  400. 

Figures  B.8-B.9  show  the  results  of  the  flow  simulation  for  the  3-D  cavity  with  3:1 
spanwise  aspect  ratio  and  Re  —  3300.  The  flow  is  computed  until  characteristic  time 
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T  —  21.  The  steady  state  has  not  yet  been  reached;  it  is  also  likely  that  a  steady 
state  may  not  exist  for  this  Reynolds  number. 

The  flow  pattern  shown  in  Fig.  B.6  in  the  symmetry  plane  is  recovered  numeri¬ 
cally  in  Fig.  B.8.  It  is  evident  that  major  recirculation  zones  are  developing  at  the 
downstream  lower  endwall,  upstream  lower  endwall  and  upper  surface  corner.  How¬ 
ever,  the  separation  at  the  upper  surface  is  not  strong  at  this  stage.  It  is  found  that 
the  viscous  boundary  layer  and  the  upper  surface  shear  layer  are  resolved  well  with 
the  clustered  (65  x  65)  grid  in  the  symmetry  plane  for  Re  =  3300.  Table  B.4  lists 
the  sizes  of  the  downstream  secondary  eddy  at  T  =  15, 17  and  21;  the  corresponding 
time-averaged  experimental  result  (Ref.  [B.  15])  is  also  included. 

Figure  B.9  shows  the  temporal  development  of  the  longitudinal  vortices  at  the 
lateral  planes  x1  =  0.5.  The  results  show  strong  3-D  effects  developing  in  the  flow. 
Two  pairs  of  TGL  vortices  have  obviously  developed  above  the  bottom  wall.  It  is 
noted  that  these  vortices  started  near  the  endwall  and  developed  towards  the  cen¬ 
tral  region.  This  indicates  that,  in  addition  to  the  downstream  concave  separation 
surface,  three-dimensionality  is  influential  in  the  formation  and  development  of  the 
longitudinal  vortex.  It  can  also  be  seen  in  these  figures  that,  to  counteract  these  TGL 
vortices,  sub- vortices  of  smaller  sizes  are  generated  very  close  to  the  wall.  Moreover, 
longitudinal  vortices  also  develop  near  the  upper  surface.  It  is  deemed  that  the 
strong  viscous  shear  layer  and  three-dimensionality  are  responsible  for  the  occurrence 
of  these  longitudinal  vortices.  It  is  observed  that  the  spanwise  grid  size  (33)  is  not 
sufficient  to  resolve  both  the  boundary  layer  near  the  endwalls  and  the  shear  layers 
forming  thereafter  in  the  interior  of  the  flow.  Evidently,  the  spanwise  grid  should  also 
be  clustered  near  the  walls. 

Examination  of  Streamwise-Vorticity  Transport 

Amongst  the  physical  phenomena  associated  with  the  3-D  flow  in  a  driven-cavity  with 
finite  span  is  the  transport  of  streamwise  vorticity,  which  is  very  important  in  helping 
understand  the  fundamental  physics  of  the  interaction  between  the  main  flow  and 
the  secondary  flows.  Therefore,  effort  was  directed  at  carefully  examining  the  flow 
structure  in  the  spanwise  direction,  in  an  attempt  to  gain  meaningful  information 
about  the  flow. 
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To  help  reduce  CPU  time  and  core  storage,  the  computational  domain  used  con¬ 
sists  of  only  one  half  of  the  cavity  (see  Fig.  B.9)  with  a  non-uniform  staggered  grid 
(65  x  65  x  49).  Symmetry  boundary  conditions  are  used  at  the  symmetry  plane  z  —  0. 
The  Reynolds  number  of  3200  is  based  on  the  velocity  (Vr)  of  the  top  surface  in  z- 
direction  and  the  depth  (Z,r)  of  the  cavity.  The  flow  is  started  impulsively  from  rest 
with  the  time-step  At  =  0.005. 

Symmetry  Plane  z  —  0 

The  2-D  flow  will  develop  in  the  symmetry  plane  z  —  0  when  a  cavity  of  in¬ 
finite  span  is  used  and  Taylor- Gortler  instability  associated  with  a  rotating  flow  is 
neglected.  Numerous  2-D  simulations  are  available.  Examination  of  the  3-D  flow 
pattern  in  this  plane  will  reveal  and  compare  the  fundamental  differences  between 
the  symmetry- plane  results  for  3-D  and  2-D  flow. 

The  results  are  illustrated  in  Fig.  B.10  using  tangential  velocity  and  contours  of 
normal  vorticity  u,  in  the  symmetry  plane,  for  the  time  instants  t  =  20,  25,  30,  35, 
45  and  50. 

As  of  t  —  20,  the  flow  pattern  in  the  symmetry  plane  shows  2-D  characteristics.  A 
primary  vortex  forms  in  the  central  region  and  three  major  secondary  eddies  near  the 
upper  back-wall,  lower  back-wall  and  lower  front-wall.  The  three-dimensional  effect 
has  not  yet  influenced  the  flow  in  this  plane. 

Between  time  t  =  20  and  t  —  25,  a  threshold  has  evidently  been  past.  The 
flow  pattern  at  t  —  25  clearly  demonstrates  the  3-D  effects.  Three  spanwise  ‘jets’  in 
the  positive  z-direction,  representing  the  Taylor-Gortler-like  (TGL)  vortices,  impinge 
onto  the  symmetry  plane.  They  interact  with  the  already- formed  primary  vorter  and 
the  three  seconday  vortices  in  the  symmetry  plane  and  drastically  change  the  flow 
pattern.  It  is  seen  that  these  three  spanswise  streams  are,  respectively,  near  the  three 
secondary  vortices  in  the  symmetry  plane,  with  the  one  near  the  lower  front-wall 
being  the  strongest.  This  suggests  that  the  three-dimensionality  is  strongest  near 
the  wall  and  three  secondary  eddies.  It  follows  that  these  effects  have  resulted  from 
the  relatively  large  shear  stresses  and  adverse  pressure  gradients  in  these  regions.  As 
the  secondary  vortex  near  the  lower  front-wall  and  the  primary  flow  above  it  are  the 
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strongest  among  the  three  secondary  eddies  in  the  symmetry  plane,  the  generation 
and  transport  of  the  ‘streamwise’  vortidty,  namely,  the  vorticity  with  its  axis  aligned 
approximately  with  the  primary  (dominant)  flow  direction,  is  expected  to  be  the  most 
significant.  The  interference  of  the  spanwise  motion  reshapes  the  flow  characteristics 
in  the  symmetry  plane,  enhancing  the  momentum  transfer  in  some  region  and  de¬ 
pressing  it  in  others. 

As  the  flow  evolves  to  t  =  30  and,  further,  to  t  =  35,  the  effects  of  three- 
dimensionality  become  more  obvious  and  significant,  particularly  near  the  lower  front- 
wall.  It  is  noted  that  the  spanwise  stream  near  the  lower  front-wall  vortex  is  moving 
towards  the  front  wall.  This  reflects  the  fact  that  the  size  of  the  secondary  vortices 
near  the  lower  front- wall  is  decreasing  with  the  increase  of  the  spanwise  motion  during 
this  period.  It  suggests  that  the  streamwise  vortices  are  generated  and  transported 
in  such  a  way  that  they  adjust  themselves  to  the  strength  of  primary  flow  and  the 
location  of  the  separation.  Experimental  results  of  Prasad  et  ai.  (1988)  have  also 
shown  this  strong  interaction  between  the  secondary  eddies  and  the  TGL  vortices. 
It  should  be  noted  that  the  mesh  in  the  central  region  of  the  symmetry  plane  is  not 
refined  enough  to  resolve  the  shear  layers  that  develop  there  at  these  two  times,  as 
shown  by  the  vorticity  contours. 

At  t  =  45,  the  3-D  pattern  is  not  as  significant  as  seen  during  t  =  25  —  35.  This 
reveals  the  unsteady  nature  of  the  TGL  vortices.  At  t  =  50,  the  3-D  flow  pattern 
takes  on  a  different  look.  The  primary  flow  near  the  bottom  wall  is  more  leveled. 
A  saddle  point  has  formed  near  the  lower  back-wall  as  a  result  of  the  interaction  of 
this  flat  primary  flow  stream,  the  back-wall  and  spanwise  motion  of  the  fluid.  The 
spanwise  motion  at  this  stage  near  the  lower  front-wall  is  not  as  strong  as  at  t  =  35. 

Although  the  above  analysis  deals  with  the  flow  structure  in  the  symmetry  plane, 
the  emphasis  has  been  stressed  on  the  interaction  of  the  spanwise  motion  (streamwise 
vortices)  with  the  flow-field  in  the  symmetry  plane.  The  flow  pattern  in  the  symmetry 
plane  is  completely  different  from  that  of  2-D  simulations.  Spanwise  motion  of  the 
fluid  evolve  and  interact  actively  with  the  primary  and  the  secondary  flows  in  the 
symmetry  plane  (as  well  as  all  x-y  planes). 

Vertical  Spanwise  Planes  x  =  0  and  x  —  .265 
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The  flow  in  the  vertical  spanwise  planes  (y-z  planes)  is  solely  due  to  three- 
dimensionality,  namely,  the  end-wall  effects  and  the  Taylor-instability  mechanism. 
Prasad  et  al.  (1988)  has  shown  that  a  sufficiently  larger  spanswise  aspect  ratio  is  also 
a  necessary  condition  for  the  existence  of  TGL  vortices.  The  experiment  of  Koseff 
and  Street  (1984)  indicated  that  unsteady  nature  of  TGL  vortices  were  obvious  for 
the  cavity  with  a  spanwise  aspect  ratio  of  3  :  1. 

For  separated  flows  with  end- walls  in  the  spanwise  direction,  it  is  known  that  "‘cor¬ 
ner  vortices”  appear  near  the  end- walls.  These  vortices  are  created  by  an  adjustment 
of  the  shear  and  pressure  forces  acting  on  the  recirculating  fluid  to  the  no-slip  condi¬ 
tion  imposed  by  the  solid  boundaries  (Ref.  [B.4]).  The  spanwise  motion  of  the  fluid 
within  the  secondary  eddies  in  the  x-y  planes  is  driven  by  these  corner  vortices.  Ob¬ 
servations  of  the  flow  pattern  in  the  symmetry  plane  in  the  previous  sub-sub-section 
have  confirmed  that  the  TGL  voriices  are  very  strong  near  the  wall  and  the  secondary 
eddies  in  the  symmetry  plane. 

Figure  B.ll  depicts  the  normal  vortidty  (a;,)  contours  in  the  vertical  spanwise 
planes  x  =  0  and  x  =  .265  for  times  t  —  20  through  t  =  50.  It  is  clearly  seen  that 
the  corner  vortices  have  originated  near  the  end- wall  z  =  -|  at  t  =  20.  The  vortices 
also  began  shedding  towards  the  center  region.  But  they  have  not  yet  penetrated 
into  the  symmetry  plane  ( z  —  0).  This  is  consistent  with  the  observation  of  the  2-D 
flow  structure  in  the  symmetry  plane  in  the  previous  sub-sub-section.  Although  some 
vortices  are  present  near  tee  symmetry  plane  at  this  stage,  they  are  not  strong  enough 
to  affect  the  flow  pattern  in  the  symmetry  plane. 

At  t  =  25,  one  and  half  pair  of  TGL  vortices  are  seen  near  the  symmetry  plane. 
Three-dimensional  effect  has  arrived  at  the  symmetry  plane.  The  2-D  pattern  in  the 
plane  has  been  disrupted.  It  is  true,  as  anticipated,  that  the  TGL  vortices  at  x  =  .265 
are  stronger  than  those  at  x  =  0  since  the  location  x  =  .265  is  in  the  vicinity  of  the 
lower  front-wall  secondary  eddy  in  the  x  —  y  plane.  It  is  also  observed  that  there  are 
some  very  weak  vortices  in  the  middle  sub-section  of  the  plane,  suggesting  that  the 
vortices  near  the  symmetry  plane  might  have  transported  from  ‘upstream’  locations 
like  x  =  .265,  where  the  TGL  vortices  developed  earlier.  It  is  interesting  to  note 
that  one  and  half  pair  of  vortices  also  appear  near  the  top  surface  and  that  their 
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axes  align  with,  those  near  the  bottom  wall  approxmately  in  the  x-y  planes  at  x  =  0 
for  t  =  25.  This  may  imply  that  the  primary  vortex  in  the  x-y  planes  possesses  a 
structure  called  a  ‘cell’  observed  experimentally  by  Mochizuki  et  al.  (1990)  at  this 
stage.  The  ‘cell’  consists  of  a  pair  of  counter-rotating  vortex  rings  of  the  primary 
vortex  and  the  norseshoe-like  secondary  vortex  ring  near  the  bottom  wall.  However, 
this  structure  is  not  obvious  at  later  times.  This  is  due  to  the  fact  that  the  experi¬ 
mental  setup  employed  a  larger  span  (aspect  ratio  as  large  as  12)  and  lower  Reynolds 
number  (Re  =  1000).  It  is  then  concluded  that  the  ‘cell’  structure  will  not  sustain 
for  the  cavity  with  an  aspect  ratio  3  :  1  at  Re  =  3200. 

The  sizes  and  the  locations  of  the  streamwise  vortices  change  with  time,  as  pre¬ 
sented  by  the  vorticity  contours  in  those  two  planes  through  t  =  50.  The  number 
of  vortices  are  increasing  in  both  planes.  The  vortices  in  the  plane  x  =  .265  seem 
stronger  than  those  in  the  plane  x  -  0.  At  the  same  characteristic  time,  there  are 
more  pairs  of  vortices  in  the  plane  x  —  .265  than  in  the  plane  x  =  0.  It  is  shown 
again  that  streamwise  vortices  near  the  secondary  eddies  in  x-y  planes  are  stronger. 
Unsteadiness  of  the  flow  at  this  Reynolds  number  Re  =  3200  is  again  vividly  illus¬ 
trated  by  the  streamwise  vorticity  changes. 

SUMMARY 

A  multigrid  distribuvive  Gauss-Seidel  (MG-DGS)  method  has  been  developed  and 
successfully  applied  to  the  simulation  of  unsteady  three-dimensional  incompressible 
flow  in  the  shear-driven  cavity.  The  MG-DGS  method  for  the  velocity  problem  has 
proved  to  be  efficient  and  robust.  The  DGS  relaxation  scheme  itself  is  very  efficient 
for  smoothing  the  high-frequency  errors  for  the  first-order  elliptic  system  for  the  ve¬ 
locity  problem.  For  uniform  grids,  convergence  rate  of  about  0.5  is  achieved  using 
accommodative  multigrid  procedures. 

Three-dimensional  shear-driven  cavity  flow  at  Reynolds  number  Re  =  3300  has 
been  simulated  using  the  presently  developed  MG-DGS  scheme  together  with  an 
ADI  algorithm.  The  longitudinal  vortices  occurring  in  cavity-flow  experiments  have 
been  captured  numerically.  Driven-cavity  flow  with  finite  span  is  definitely  three- 
dimensional.  End-walls  play  a  dominant  role  in  creating  streamwise  vortices  (TGL 
vortices),  The  streamwise  TGL  vortices  interact  strongly  with  primary  and  secondary 
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vortices  in  the  x-y  planes.  The  evolution  and  transport  of  streamwise  vorticity  alters 
the  structure  and  sizes  of  the  primary  and  secondary  vortices  in  the  symmetry  plane 
(and  all  other  x-y  planes).  The  flow  is  persistently  unsteady  for  Re  =  3200  through 
t  —  50.  The  flow  structure  observed  in  the  numerical  simulation  is  in  good  agreement 
with  available  experimental  observations. 
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2.C 


Adaptive  Grid  Technique  for  Higher-Re  Flows 


2.C.1  Abstract 

The  objective  of  this  study  is  to  efficiently  simulate  vortex-dominated  highly  unsteady  flows. 
In  such  flows,  the  locations  as  well  as  the  extent  of  the  regions  requiring  fine-mesh  resolution  vary 
with  time.  A  technique  has  been  developed  to  simulate  these  flows  on  a  temporally  adapting  grid 
in  which  the  adaption  is  based  on  the  evolving  flow  solution.  The  flow  in  an  axisymmetric 
constriction  has  been  selected  as  an  illustrative  problem.  The  multiple  and  disparate  length  scales 
inherent  in  this  complex  flow  make  this  problem  ideally  suited  for  evaluating  the  adaptive-grid 
technique.  Adaption  is  based  on  the  equidistribution  of  a  weight  function,  through  the  use  of  forcing 
functions.  The  significance  of  this  is  that  the  method  can  be  implemented  into  existing  flow-analysis 
systems  with  minimal  changes.  The  grid-generation  equations  developed  are  viewed  as  grid-transport 
equations.  The  time-dependent  control  functions  perform  the  role  of  the  convective  speed  in  this 
transport  mechanism.  The  equations  provide  the  efficiency  and  flow  tracking  capability  of  parabolic 
equations,  while  maintaining  the  smoothness  of  computationally  expensive  elliptic  equations.  The 
efficiency  and  flow  tracking  capability  of  the  approach  is  demonstrated  for  both  steady  and  unsteady 
flows. 

2.C.2  Background 

The  use  of  appropriate  coordinates  in  the  mathematical  formulation  of  a  problem  has  a 
dominant  effect  on  the  accuracy  and  computational  resource  requirements  of  analysis  and  numerical 
simulation  of  fluid  flow  problems.  In  addition  to  requiring  that  the  coordinates  be  boundary- 
oriented,  it  is  desirable  that  they  also  reflect  some  of  the  significant  features  of  the  flow.  Traditional 
grid-generation  methods  rely  on  the  experience  of  the  user  and  only  partial  a  priori  knowledge  of 
localized  or  transient  physical  phenomena  which  require  a  fine  mesh  for  their  accurate  representation. 
The  objective  of  the  present  work  is  to  develop  an  efficient  flow-simulation  technique,  which 
automatically  resolves  critical  regions  of  widely  disparate  length  scales.  A  major  advantage  of  a 
time-dependent  mapping  incorporating  the  evolving  flow  physics  is  its  ability  to  cluster  coordinate 
lines  only  where  necessary  at  that  instant,  thus  reducing  meshpoint  requirements.  This  feature  is 
particularly  useful  for  strongly  unsteady  flows,  in  which  the  same  flow  regions  do  not  require  fine 
resolution  at  all  instants  of  time.  For  flows  requiring  fine  resolution  throughout  the  domain,  such 
as  intensely  mixing  turbulent  flows,  the  advantage  is  greatly  diminished,  if  not  eliminated.  Several 


review  papers  such  as  by  Thompson  (1985),  Eiseman  (1985)  and  Anderson  (1987)  present 
applications  where  significant  improvements  in  accuracy  and  efficiency  have  been  obtained  through 
the  use  of  adaptive-grid  procedures.  Many  of  the  procedures  require  the  solution  of  elliptic  systems 
for  each  grid  regeneration,  as  done,  for  example,  by  Kim  and  Thompson  (1990)  and  in  the 
variational  approach  of  Brackbill  and  Saltzman  (1982),  rendering  the  procedure  rather  costly  for  truly 
transient  problems.  Bockelie  (1988)  and  Trompert  (1990)  have  developed  techniques  for  transient 
problems  based  on  iterative  procedures  requiring  multiple  grid  subdivisions  for  each  grid 
regeneration  or  time  level  advancement.  The  present  adaptive-grid  technique  draws  on  the  error- 
equidistribution  concept  of  Anderson  (1987),  together  with  the  widely  used  Poisson  generation 
system  of  Thompson  (1974),  with  the  important  modification  that  the  equations  governing  the  grid 
evolution  are  parabolic  in  time.  These  grid-transport  equations  fit  well  into  the  flow  solution 
procedure  and  allow  grid  regeneration  at  every  time  level,  with  reasonable  computational 
requirements. 

In  the  following  sections,  the  equations  governing  the  flow  are  expressed  in  generalized  time- 
dependent  coordinates,  and  the  grid-transport  equations  and  weight  functions  are  derived.  The 
numerical  procedure  is  described  briefly  and  results  are  presented  and  discussed  for  a  model  problem 
of  an  axisymmetric  internal  flow. 

2.C.3  Flow  Governing  Equations  in  Generalized  Time-Dependent 
Coordinates 

The  present  analysis  considers  time-dependent  flows  of  incompressible,  Newtonian  and 
constant-viscosity  (v)  fluids.  Therefore,  the  incompressible  unsteady  Navier-Smkes  equations  are 
the  appropriate  governing  equations.  Written  in  terms  of  the  vorticity  vector  u  and  the  velocity 
vector  v ,  these  equations  consist  of  a  temporally  parabolic  vorticity-transport  equation  given  as 

+  (v-V)<3  =  ( (3 ’V)  v  ~  (VxVxcS)  /  (C.l) 

dt  Re 


with  the  following  kinematic  definition  for  vorticity: 

(3  =  Vxf?.  (C.2) 

Equations  (C.l)  and  (C.2)  have  been  non-dimensionalized  using  LR  as  the  characteristic  length,  UR 
as  the  characteristic  speed,  LR/UR  as  the  characteristic  time,  and  UR/LR  as  the  characteristic  unit  of 
vorticity.  The  Reynolds  number.  Re,  is  based  on  UR  and  LR  to  be  defined  later  for  the  flow  problem 
studied. 


The  axisymmetric  flows  considered  involve  only  two  independent  spatial  dimensions.  Hence, 
the  mass-conservation  equation  is  identically  satisfied  by  the  introduction  of  the  Stokes  stream 
function,  i p,  related  to  the  local  velocity  vector  as 


v  =  VjJt  x  e3 , 


(C.3) 


where  e 3  is  the  fundamental  contravariant  base  vector  locally  normal  to  the  (x,r)  meridional  plane, 
(Fig.  1). 

The  governing  equations  (C.1)-(C.3)  are  in  general  vector  form.  For  application  purposes, 
a  coordinate  system  must  be  selected  and  Eqs.  (C.1)-(C.3)  represented  in  scalar  component  form. 
The  boundary-oriented  time-dependent  ‘computational  coordinates’  ($l,?2,f\r)  to  be  employed  are 
defined  through  an  admissible  coordinate  transformation  T.  The  implementation  of  axisymmetry  is 
simplified  by  using  a  transformation  of  the  following  form, 

x1  =  x 1  (E1  l2  t ) 

x2  =  zal.V.*)co3[V)  (C.4) 

x3  =  r  (l1,  l2.  x)  sin($3)  . 


The  physical  space  is  represented  by  the  basic  Cartesian  coordinate  system,  (x',x2,xJ)  or  (x,y,z), 
while  and  £2  represent  the  generalized  curvilinear  coordinates  in  a  given  meridional  plane,  and  £3 
represents  the  transverse  or  azimuthal  coordinate  <t>  describing  the  various  meridional  planes,  (Fig. 
C.l). 


In  terms  of  the  stream  function  defined  in  Eq.  (C.3),  the  definition  of  vorticity  given  by  Eq. 
(C.l)  expressed  in  generalized  coordinates  (Osswald  [1983])  becomes: 


a_ 
Jg'K1 


-it 


sfg  dV  Sg  as1 


=  -Q3 


(C*  5) 


where 


g  =  det  . 


(C.6) 


The  defining  equation  (C.5)  for  w3  constitutes  an  elliptic  partial  differential  equation  governing  the 
development  of  the  axisymmetric-flow  stream  function. 


The  vorticity-transport  equation,  Eq.  (C.  1),  expressed  in  terms  of  the  generalized  coordinates 
and  written  in  strong  conservation-law  form,  is  as  follows  (Thornburg  [1991]): 

Equations  (C.5)-(C.7)  comprise  the  flow  governing  equations  in  generalized  time-dependent 
coordinates.  The  next  section  presents  the  grid-transport  equations  used  to  define  the  coordinate 
transformation  given  by  Eq.  (C.4). 


2.C.4  Grid-Transport  Equations 


To  facilitate  the  determination  of  accurate  flow  solutions,  the  physical  domain  of  arbitrary 
shape  is  transformed  into  a  computational  domain  of  regular  rectangular  shape,  such  that  a  uniform 
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(C.7) 


computational  mesh  corresponds  to  an  appropriately  distributed  physical  grid.  The  goal  is  to 
incorporate  the  developing  solution  variables  into  the  evolving  transformation  such  that  the  available 
grid  points  are  distributed  in  a  manner  that  represents  the  physical  solution  with  the  ‘best’  possible 
accuracy.  This  leads  to  the  search  for  the  relationship  between  the  developing  physical  solution  and 
the  corresponding  grid  distribution.  Unlike  the  well  defined  conservation  laws  governing  fluid  flow, 
there  is  no  unique  physical  law  governing  the  relationship  between  the  solution  and  the  best  grid. 

The  requirement  of  strong  smoothing  properties  and  tendency  towards  uniformity  suggests 
that  a  law  elliptic  in  nature  would  be  most  appropriate.  However,  Laplace’s  equation,  the  simplest 
elliptic  equation,  provides  no  control  over  grid  spacing.  Grids  generated  by  Poisson  equations 
achieve  non-uniform  spacing  enforced  by  the  non-homogeneous  terms,  called  forcing  functions  in 
this  context;  for  examples,  see  Thompson  (1985b).  To  develop  a  time-dependent  transformation, 
these  forcing  functions  can  be  based  on  the  evolving  flow  solution.  For  each  time  step,  grid 
regeneration  would  require  the  solution  of  an  elliptic  system.  This  does  not  seem  practical  for  a  truly 
unsteady  flow  situation. 

For  unsteady  flow  problems,  such  as  the  one  considered  in  the  present  study,  the  flow- 
transport  equations  are  parabolic  in  time  and  elliptic  in  space.  Thus,  a  parabolic  law  governing  the 
grid  motion,  that  could  be  integrated  in  time  along  with  the  flow  transport  equations,  would  fit  well 
with  the  solution  procedure.  Grid  generation  methods  based  on  parabolic  systems  provide 
computational  efficiency,  but  often  sacrifice  grid  smoothness  and  control  over  clustering. 

Ghia,  Ghia  and  Shin  (1983)  developed  a  grid  generation  system  by  setting  to  zero  the 
coefficient  of  the  convective  term  in  the  two-dimensional  Navier-Stokes  equations  in  the  transformed 
plane.  The  equations  thus  obtained  have  the  form  of  a  transport  equation.  Their  procedure  then 
employed  the  time-dependent  method  in  order  to  march  both  the  flow  and  the  grid  equations  to 
steady  state.  In  that  work,  the  flow  velocities  were  used  to  drive  the  grid  movement.  It  has  been 
recognized  (Anderson  [1987])  that  the  Poisson  equation  of  Thompson  (1974)  in  the  transformed 
plane  can  be  rearranged  into  a  transport-equation  form  by  the  addition  of  a  temporal  derivative. 
Therefore,  the  grid  equations  are  written  as  follows: 

The  time-dependent  control  functions  perform  the  role  of  the  convective  speed  for  the  grid-transport 
mechanism.  This  allows  the  grid  to  be  advanced  in  time,  along  with  the  parabolic  flow-transport 
equation.  The  grid  equations  are  dynamically  coupled  with  the  evolving  flow  field  through  the  P  and 
Q  terms. 
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J2  ( Px{  +  Qx^  +  xj 
J2[Pr{  +  C?rn  +  rj 


-  2Px?n  +  y^Sv 
ar«  ~  2P^n  +  Yr,,. 


(C.8) 


Approximate  Equidistribution  through  Forcing  Functions 

The  equidistribution  of  a  weight  function  on  a  mesh  is  the  main  concept  of  the  present 
adaptive  grid  scheme.  It  has  been  recognized  (Anderson  [1987])  that  a  one-dimensional 
equidistribution  law  may  be  written  as  a  Poisson  system, 

Xtt  *  Px f  =  0,  with  P  -  .  (C .  9 ) 

w 

Using  variational  calculus,  equidistribution  laws  may  be  similarly  written  for  higher  dimensions.  In 
order  to  simplify  the  equations  and  reduce  the  computational  requirements,  the  problem  is  viewed 
as  a  series  of  one-dimensional  equidistribution  schemes.  The  individual  forcing  functions  are 
combined  to  yield  the  forcing  function  for  multiple  dimensions  as 


Equations  (C.8)  and  (C.10)  provide  an  approximate  equidistribution  law  based  upon  the  individual 
weight  functions.  The  development  of  the  weight  functions  for  the  individual  coordinate  directions 
is  described  in  the  following  section. 

Derivation  of  Weight  Functions 

Grid  spacing  is  inversely  proportional  to  the  weight  function,  and  hence,  the  weight  function 
determines  the  grid  distribution.  The  weight  function  used  in  the  present  work  is  an  approximation 
to  the  truncation  error.  Determination  of  this  function  is  one  of  the  most  challenging  areas  of 
adaptive  grid  generation.  It  is  believed  that  the  function  developed  in  the  present  work  represents 
a  significant  contribution  in  this  area.  For  a  k*-order  accurate  simulation,  the  leading  term  in  the 
truncation  error  for  convective-like  derivatives  is  of  the  form  hk  0lk+!),  where  h  is  the  local  mesh 
spacing,  0  is  the  solution  variable  and  0(k+I)  is  the  (k+ 1)*  derivative  of  0;  similarly  for  diffusion-like 
derivatives  it  is  of  the  form  hk  0(k+2>.  If  a  flow  simulation  is  to  benefit  from  non-uniform  'pacing, 
the  flow  must  contain  regions  of  non-uniform  derivatives.  An  equidistributional  scheme  will  result 
in  non-uniform  spacing  when  0(k+1)  varies  over  the  solution  domain.  Non-uniform  spacing  will  lead 
to  the  greatest  benefit  when  the  variation  in  0(k+ 1)  is  largest.  This  suggests  the  use  of  derivatives 
in  developing  the  weight  function.  Evaluation  of  higher-order  derivatives  from  discrete  data  is 
progressively  less  accurate  and  subject  to  noise.  However,  variation  in  the  derivative  0(k'H)  implies 
that  lower-order  derivatives  also  vary.  In  fact,  if  0Ck+1)  varies  rapidly,  0k  can  be  large.  When  0(k+1) 
varies  widely,  the  lower-order  derivatives  must  be  non-zero  in  the  vicinity  of  wide  variation,  and  are 
proportional  to  the  rate  of  variation.  Therefore,  it  is  possible  to  employ  some  lower-order 
derivatives  as  a  proxy  for  the  truncation  error. 
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Experience  with  flow  simulation  also  sheds  light  on  flow  regions  requiring  fine  resolution. 
Thus,  Flow  regions  of  large  gradients  require  finer  resolution  than  more  uniform  flow  regions,  as 
can  be  deduced  from  the  above  discussion.  Also,  regions  of  large  solution  variable  curvature,  which 
necessarily  occur  between  regions  of  high  and  low  gradients,  require  resolution.  Many  researchers 
have  successfully  employed  weight  functions  composed  of  solution  variable  gradient  and  curvature 
(Anderson  [1987],  Thompson  [1985]).  Ghia,  Ghia  and  Shin  (1983)  have  clearly  demonstrated  the 
importance  of  properly  resolving  regions  of  large  solution  variable  curvature  for  some  model 
problems.  The  present  weight  functions  employ  higher-order  derivatives  of  the  physical  velocity 
variables  u  and  v  through  the  derived  variable  of  vorticity.  The  choice  of  vorticity  is  logical,  as 
vorticity  provides  a  physical  description  of  the  rotational  intensity  of  the  flow.  One  difficulty  in  the 
use  of  vorticity  for  constructing  the  weight  function  is  that  it  may  vary  over  several  orders  of 
magnitude.  This  can  result  in  a  grid  equation  system  that  is  very  stiff  and,  hence,  difficult  to  solve 
numerically  and,  in  fact,  has  historically  limited  its  application.  In  addition,  important  information 
can  be  lost  as  derivatives  in  regions  of  large  vorticity  swamp  important  derivatives  in  regions  of 
lower  vorticity,  much  as  the  ‘convective’  terms  of  a  transport  equation  can  swamp  out  important 
diffusion  terms  at  higher  Reynolds  number.  The  newly  derived  weight  functions  overcome  these 
difficulties  by  an  appropriate  normalization  procedure,  resulting  in  well  behaved  grid  equations  as 
well  as  adaption  in  lower  vorticity  regions  that  demand  fine-mesh  resolution  for  accurate  numerical 
simulation. 


The  weight  functions  have  been  developed  using  a  component-by-component  approach.  Each 
component  has  been  designed  to  increase  grid  clustering  in  regions  where  a  certain  flow  situation 
occurs.  The  component-by-component  approach  also  facilitates  the  inclusion  of  additional  terms,  if 
the  need  arises  in  future  studies.  The  weight  functions  formulated  are  as  follows: 
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The  first  term  is  the  constant  1.0  and  is  employed  to  ensure  that  the  weight  function  is  always  non¬ 
zero  positive.  A  zero  weight  function  would  imply  infinite  grid  spacing.  If  the  weight  function  were 
to  become  negative,  equidistribution  would  require  a  negative  volume  of  the  computational  cell.  The 
second  term  is  the  normalized  magnitude  of  the  vorticity  and  results  in  grid  clustering  in  regions  of 
large  vorticity,  or  large  velocity  gradients.  This  term  does  not  account  for  resolution  needed  in 
regions  of  large  vorticity  gradient  or  slope  but  small  vorticity  magnitude.  Such  a  region  can  be  a 
source  of  significant  truncation  error.  In  fact,  this  is  often  the  case  in  unsteady  or  higher  Reynolds 
number  flows. 


The  second  term  is  the  normalized  relative  derivative  of  the  vorticity  with  respect  to  the 
coordinate  £k.  The  contribution  of  this  term  is  from  regions  of  large  vorticity  gradients  that  need 


49 


resolution,  such  as  those  adjoining  high-vorticity  regions.. 

The  third  term  is  the  normalized  relative  second-order  derivative  of  vorticity  with  respect  to 
the  coordinate  £k.  The  use  of  this  term  is  an  attempt  to  resolve  regions  of  large  vorticity  curvature 
near  the  periphery  of  large  vorticity-gradient  regions.  The  use  of  the  relative  derivative  is  important 
to  eliminate  domination  by  large  vorticity  gradients  present  in  high-vorticity  regions.  Thus,  the 
relative  derivative  is  a  better  measure  of  the  local  truncation  error  than  the  absolute  derivative. 

Evaluation  of  Weight  Function  Coefficients 


The  weight  functions  are  determined  by  the  vorticity  distribution  as  well  as  by  the  coefficients 
ak,  bk,  and  ct,  in  Eq.  (C.  11).  A  wide  range  of  values  for  these  coefficients  has  been  tested.  The 
choice  of  values  has  been  based  primarily  on  numerical  experimentation.  In  general,  a  fairly  wide 
range  of  values  for  these  parameters  yields  acceptable  results. 

The  contribution  of  each  weight-function  component  is  determined  by  its  coefficients.  A  10 
percent  contribution  from  the  magnitude  of  the  vorticity,  50  percent  from  the  relative  first-derivative 
and  40  percent  from  the  relative  second-derivative  of  vorticity  has  been  found  to  provide  the  best 
resolution  for  the  flows  studied.  This  distribution  determines  the  relative  magnitude  of  a1,  a2,  bl, 
b2,  c1,  and  c2.  Theoretically,  the  ratio  of  the  minimum  value  of  the  weight  function  to  its  maximum 
value  is  equal  to  the  ratio  of  maximum  cell  volume  to  minimum  cell  volume.  The  present  method 
is  approximately  equidistributional  and,  as  such,  this  relation  is  not  exact.  However,  estimates  of 
the  maximum  and  minimum  values  of  the  weight  functions  can  be  used  to  predict  the  ratio  of 
maximum  to  minimum  cell  volume.  Hence,  a  limit  may  be  placed  on  the  maximum  allowable  cell 
size  variation.  The  present  work  uses  a  value  of  200,  so  that 
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The  massive  grid  migration  necessary  to  maintain  this  theoretically  optimum  grid  distribution  in  the 
streamwise  direction  requires  a  time  step  much  smaller  than  is  needed  to  resolve  the  temporal  scales 
of  the  physical  problem.  This  is  not  efficient.  A  trade-off  exists  between  larger  values  and  greater 
adaption,  versus  overall  efficiency.  Stronger  adaption  produces  greater  grid-point  migration, 
resulting  in  slower  convergence  of  the  stream- function  equation,  as  the  initial  guess  is  then  further 
away  from  the  true  solution.  A  reduction  in  time  step  may  also  be  required  for  very  large  grid 
movement.  Thus,  with  the  present  adaptive  method,  and  likely  with  others,  the  smallest  possible 
number  of  grid  points  may  not  necessarily  lead  to  the  greatest  efficiency.  For  flows  with  widely 
spaced  and  rapidly  moving  flow  structures  requiring  resolution,  such  as  a  vortex-ring  convecting 
downstream,  it  is  more  efficient  to  use  a  few  additional  grid  points  and  moderate  adaption.  Thus, 
the  parameters  should  be  chosen  such  that  the  time  step  limitation  is  imposed  by  the  time  scales  of 
the  flow,  and  not  by  the  degree  of  grid  movement.  These  comments  apply  primarily  to  the 
streamwise  direction,  where  the  large  extent  of  the  domain  permits  grid  points  to  travel  large 
distances.  The  confined  nature  of  the  geometry  in  the  normal  direction,  combined  with  the  spatially 
elliptic  equations  governing  grid  movement,  tends  to  moderate  movement  in  the  normal  direction. 
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Hence,  for  the  current  work,  the  coefficients  are  scaled  by  a  factor  of  1000  in  the  streamwise 
direction.  The  values  employed  for  the  weight  function  coefficients  are  a1  =  20,  b1  =  100,  c1  = 
80,  a2  =  .02,  b2  =.l  and  c2  =.08.  For  the  family  of  geometries  studied,  simulations  indicate  that 
the  values  of  these  coefficients  can  be  held  constant,  independent  of  Re,  Ar  or  streamwise  extent  of 
physical  domain. 

The  most  important  coefficient  in  the  expression  for  the  weight  function  is  that  of  the  first- 
derivative.  The  contribution  of  the  vorticity  magnitude  term  is  primarily  in  regions  of  large  flow 
acceleration  and/or  shear,  such  as  near  the  throat.  The  importance  of  this  term  decreases  as  the 
unsteadiness  of  the  flow  increases.  The  first-derivative  term,  in  conjunction  with  the  magnitude,  has 
been  sufficient  to  resolve  all  but  the  most  severe  cases.  For  the  latter,  additional  resolution  was 
required  in  regions  of  large  vorticity  curvature.  Therefore,  the  second-order  derivative  term  was 
implemented  to  successfully  resolve  these  regions.  This  is  particularly  important  for  the  near-wall 
region,  where  resolution  is  important  but  both  magnitude  and  first  derivatives  of  the  solution  flow 
variables  may  be  moderate.  Also,  this  term  was  found  to  be  crucial  for  appropriate  resolution  of  the 
region  of  low-speed  flow  along  the  wall  and  upstream  of  the  constriction  and,  to  avoid  contamination 
of  the  flow  downstream  of  this  region. 

2.C.5  Model  Internal  Flow  Problem 

The  geometry  of  the  test  problem  chosen  was  a  family  of  axisymmetric  constrictions 
described  by  a  cosine  function,  and  is  illustrated  schematically  in  Fig.  C.2.  For  this  geometry, 
experimental  data  is  available  for  the  Reynolds  number  range  500-15,000  and  for  constrictions  of 
25-75  percent  by  area  (Ahmed  and  Giddens  [1983]  and  Deshpande  and  Giddens  [1980]).  The 
asymptotic  state  of  the  flow  is  steady  or  unsteady,  depending  upon  the  severity  of  the  constriction 
and  the  value  of  the  Reynolds  number. 

Numerical  Solution 

Equations  (C.5),  (C.7)  and  (C.8)  are  written  in  computational  coordinates.  The  implied 
transformation  maps  the  general  region  of  interest  R  in  the  physical  (x,r)  plane,  to  the  unit  square 

□  in  the  computational  plane  defined  as  n  a  [(£,??)  |  0<  £<1,  0<tj^1].  A  uniformly  spaced 
finite  difference  grid  A  of  size  (N  x  M)  is  employed  to  discretize  the  unit  computational  region 

□  . 


The  three  coupled  equations,  two  for  the  grid  and  one  for  the  vorticity,  are  integrated  in  time. 
The  grid  equations  share  some  of  the  ‘stiffness’  of  the  flow  transport  equations  through  the  velocity¬ 
like  forcing  functions  P  and  Q.  In  fact,  the  transformation  may  be  viewed  as  transferring  the 
difficulty  usually  encountered  in  solving  the  flow  equations,  to  that  now  experienced  in  solving  the 
grid  equations.  The  ‘stiffness’  of  the  flow  equations  is  reduced  through  the  scaling  of  the 
transformation  metrics  appearing  in  their  coefficients.  The  equations  are  discretized  and  solved 
numerically.  Time  derivatives  are  approximated  using  a  first-order  accurate  backward  difference. 
The  discretized  non-linear  terms  are  quasi-linearized.  Extrapolation  in  time  is  employed  to  generate 
an  initial  guess  for  the  iterative  solution  of  the  stream  function.  In  addition,  a  local  iteration 
procedure  is  employed  to  update  the  coefficients  of  the  ‘convective’  terms  in  the  vorticity-transport 
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equation.  For  further  details  of  the  procedure,  see  Thornburg  (1991).  The  location  and  spacing  of 
the  mesh  points,  not  the  accuracy  with  which  the  grid  equations  are  evaluated,  affects  the  accuracy 
and  stability  of  the  discretized  flow  equations.  Therefore,  in  order  to  produce  a  numerically  stable 
solution  to  the  grid  equations,  the  idea  of  upwind  differencing  is  used  in  these  equations.  Upwind 
differencing  of  the  flow  equations  adversely  affects  the  accuracy  of  the  flow  solution  through  the 
introduction  of  numerical  diffusion,  which  lowers  the  local  effective  Reynolds  number.  However, 
this  very  property  is  beneficial  for  the  grid  equations.  It  is  desirable  for  the  grid  to  vary  smoothly 
near  regions  of  large  gradients  where  grid  clustering  occurs.  Thus,  the  numerical  diffusion 
introduced  by  first-order  upwind  differencing  of  the  grid  equations  not  only  adds  to  the  stability  of 
the  grid-transport  process,  but  actually  improves  the  grid  quality,  and,  hence,  the  flow  simulation 
accuracy.  The  biased  differencing  is  based  on  the  velocity-like  P  and  Q  forcing  functions. 


The  source  terms  of  the  discretized  equations  require  the  values  of  P  and  Q  at  discrete 
locations,  and  must  aiso  be  evaluated  from  discrete  weight  function  data.  The  weight  functions 
contain  values  of  the  vorticity  field  at  time  level  (n+1)  and  represent  the  dynamic  coupling  between 
the  grid  and  the  flow  solution.  The  simplest  technique  to  decouple  the  equations  would  be  to  lag  the 
vorticity  appearing  in  the  grid  forcing  functions.  This  approach  was  implemented  and  found  to  be 
less  than  optimal.  A  very  small  time  step  was  needed  to  adequately  track  the  temporally  evolving 
flow  even  for  moderate  Reynolds  number  cases.  Another  contribution  of  the  present  work  lies  in 
addressing  this  shortcoming,  by  using  temporal  extrapolation  to  obtain  a  guessed  value  for  vorticity 
as: 


o>T.j  =  2o )nitj  -  «?;}  +  O  ( At) 2 .  (C.  13 ) 

These  extrapolated  values  for  the  vorticity  are  used  to  evaluate  the  forcing  functions. 

High-Reynolds  number  flows  often  exhibit  localized  regions  of  intense  variation  of  vorticity. 
Hence,  the  variation  of  the  weight  function  can  become  severe  and  localized  in  this  region,  resulting 
in  rapid  localized  variation  in  grid  spacing.  The  interval  of  influence  of  these  localized  regions  is 
increased  by  smoothing  the  forcing  functions,  resulting  in  smoother  variation  of  mesh  size. 

In  order  to  control  the  truncation  error  associated  with  the  grid  speed  terms,  as  well  as 
provide  stability  for  the  grid  evolution,  it  is  necessary  to  control  the  magnitude  of  grid  movement. 
The  unsteady  nature  of  the  present  flow  at  higher  Re  produces  vortical  structures  which  convect  out 
of  a  finite  solution  domain.  This  can  cause  excessive  and  unnecessary  grid  movement  at  the  outlet, 
because  as  a  structure  leaves  the  solution  domain,  its  associated  grid  points  rapidly  migrate  to  the 
next  structure  upstream.  In  practice,  it  was  found  beneficial  to  create  a  buffer  region  at  the  outlet 
for  the  grid  solution.  The  buffer  region  is  created  by  exponentially  damping  the  forcing  functions 
to  zero  over  the  last  25  points  of  the  domain  in  the  streamwise  direction.  The  buffer  region  is  used 
to  eliminate  the  rapid  grid  movement  in  this  region  of  relatively  quiescent  flow,  thus  removing  the 
time  step  constraint  imposed  by  this  movement. 

Overhead  Associated  with  Adaption 
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For  the  unsteady  flow  cases  examined,  the  adaption  procedure  is  enabled  throughout  the 
simulation.  By  monitoring  the  residence  time  of  each  subroutine  associated  with  the  adaption,  the 
direct  overhead  is  found  to  be  of  the  order  of  15-30  percent  for  the  unsteady  cases  studied.  The 
overhead  for  steady- flow  cases  is  much  less,  as  the  adaption  procedure  is  active  for  only  a  small 
fraction  of  the  total  simulation  time.  The  iterative  method  employed  for  the  steady  problems  studied 
exhibits  rapid  convergence  initially.  The  grid-transport  equations  do  not  represent  physical  laws,  and 
hence,  the  level  of  convergence  required  for  acceptable  accuracy  is  much  less  than  that  for  the  flow- 
transport  equations.  For  steady-flow  or  even  quasi-steady  flow,  no  further  accuracy  is  achieved  in 
the  flow  solution  by  adaption  after  the  general  flow  structure  has  developed  and  the  grid  adapted  to 
it.  Hence,  for  steady-flow  calculations,  the  grid  movement  is  stopped  before  a  fully  converged  flow 
solution  has  been  achieved.  For  a  steady  problem  where  the  initial  transients  are  not  of  interest, 
it  seems  wasteful  to  spend  resources  to  adapt  and  resolve  these  transients.  Therefore,  the  most 
efficient  procedure  for  steady-flow  problems  has  been  found  to  consist  of  starting  the  calculation  with 
some  initial  static  grid  and  running  the  simulation  until  the  gross  structure  has  developed.  Then  the 
adaption  procedure  is  activated  for  approximately  50-100  iterations,  and  deactivated  thereafter  as  grid 
movement  becomes  small.  This  is  also  a  sufficient  number  of  time  integrations  for  the  adaption 
procedure  to  sense  and  automatically  resolve  flow  structures  appearing  due  to  the  increased  resolution 
provided  by  the  adapting  grid  in  certain  regions.  The  simulation  is  then  continued,  using  the  static 
but  adapted  grid,  to  the  desired  level  of  convergence  for  the  flow  solution.  This  approach  greatly 
reduces  the  overhead  of  the  adaptive  procedure  for  steady-flow  problems. 


Resolution  of  Scales 

The  resolution  of  the  length  and  time  scales  prevailing  in  a  flow  is  the  most  important  factor 
for  an  accurate  simulation.  Time-accurate  simulations  of  unsteady  flow,  such  as  those  performed 
in  the  present  work,  must  employ  a  spatially  constant  time  step.  For  the  flows  studied,  it  is  believed 
that  the  reduction  in  the  required  number  of  grid  points  due  to  adaption  is  greatest  in  the  streamwise 
direction,  where  the  reduction  has  been  as  much  as  50  percent.  The  reduction  of  points  in  the 
normal  direction  is  not  as  dramatic.  The  regions  of  intense  flow  gradients  requiring  normal 
resolution  travel  over  a  shorter  distance  and  their  locations  are  easily  predicted.  For  example,  the 
oscillating  shear  layer,  which  requires  the  greatest  normal  resolution,  is  confined  to  a  location  only 
several  times  its  thickness.  In  addition,  the  wall  region  always  requires  normal  resolution.  Using 
this  a  priori  information,  an  initial  grid  with  a  fairly  efficient  normal  grid-point  distribution  can  be 
constructed,  lessening  the  need  for  an  adaptive  procedure  in  this  coordinate  direction.  Thus,  the 
reduction  in  the  number  of  grid  points  required  in  the  normal  direction  is  estimated  to  be  20-30 
percent,  as  compared  with  a  static  clustered  grid,  depending  on  the  particular  case. 

Proper  resolution  of  the  time  scales  was  a  much  more  difficult  problem  than  anticipated.  The 
time  scales  of  the  flows  in  the  more  severe  geometries  were  very  short.  For  the  configurations 
studied,  the  severity  of  the  flow  was  primarily  due  to  Uie  geometric  configuration.  The  most 
important  consideration  being  the  area  constriction  ratio  (ACR),  followed  by  the  half-length  (XSL), 
over  which  the  constriction  occurred.  Plots  of  the  truncation  error  associated  with  the  time 
derivative,  based  on  static  grid  cases,  showed  that  the  shear  layer  immediately  downstream  of  the 
obstruction  exhibited  the  largest  of  these  truncation  errors  and  the  shortest  time  scale. 
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It  is  not  straightforward  to  analyze  the  truncation  errors  associated  with  time  derivatives  for 
moving-  grid  simulations  and,  hence,  to  determine  if  the  At  used  provides  sufficient  temporal 
resolution  for  the  method  employed.  The  most  reasonable  method  for  determining  temporal 
resolution  consists  of  carrying  out  several  simulations  with  various  time  steps,  constant  for  a  given 
simulation.  Time-step  independence  suggests  sufficient  temporal  resolution.  For  ACR  =  0.75  with 
XSL  =  2.0,  the  flow  with  Re  =  2000  has  been  analyzed  extensively.  This  case  was  chosen  for 
study  due  to  the  fact  that  much  experimental  data  was  available  for  it.  A  simulation  using  first-order 
accurate  time  discretization  was  performed,  and  progressively  smaller  time  steps  were  employed  until 
successive  calculations  showed  no  qualitative  variation.  Time  step  independence  was  observed  for 
At  <  0.001.  In  addition,  a  simulation  using  second-order  accurate  time  differencing  was  also 
performed  for  this  case.  No  qualitative  difference  was  observed  between  the  two  sets  of  results. 
It  is  therefore  believed  that  the  first-order  accurate  temporal  discretization  with  a  time  step  of  0.001 
provides  adequate  temporal  resolution  for  this  flow  situation.  This  procedure  is  similar  to  varying 
grid  sizes  to  evaluate  grid  independence  and,  hence,  spatial  resolution.  For  less  severe  geometries 
such  as  the  one  with  ARC=  0.50,  XSL  =  2.0  and  Re  =  2000,  the  time  scales  are  much  less  severe. 
For  this  case,  no  qualitative  differences  were  observed  in  the  results  using  either  first-order  accurate 
or  second-order  accurate  temporal  discretization,  so  long  as  the  time  step  employed  was  smaller  than 
0.002.  However,  small  scale  differences  were  observed  under  quantitative  analysis. 

2.C.6  Simulation  Results  and  Discussion 

Results  are  presented  for  various  constriction  ratios  and  Reynolds  numbers.  The  severity  of 
the  geometry  triggers  separation  at  very  low  Reynolds  number.  For  example,  for  the  configuration 
with  ACR  =  0.75  and  XSL  =  1.0,  separation  occurs  at  a  Reynolds  number  of  10  (Lee  and  Fung 
[19701).  The  flows  simulated  are  grouped  in  two  categories.  The  first  category  deals  with  steady 
flows,  while  the  second  deals  with  unsteady  flows.  The  category  to  which  a  given  flow  belongs 
depends  on  the  values  of  the  three  parameters,  namely,  area  constriction  ratio  (ACR),  physical  half- 
length  (XSL)  of  constriction  and  Reynolds  number  (Re). 


Steady-Flow  Results 


A  series  of  calculations  were  performed  in  order  to  validate  solutions  on  non-uniform,  non- 
orthogonal  grids  and  to  demonstrate  grid  independence.  The  geometrical  configuration  with  ACR 
=  0.50,  XSL  =  2.0  and  Re  =  500  was  chosen,  due  to  the  availability  of  experimental  data  for  it. 
The  first  result  presented  is  a  simulation  with  a  large  (450  x  45),  essentially  uniform,  orthogonal 
grid,  and  was  undertaken  to  obtain  a  benchmark  set  of  data.  The  grid  as  well  as  the  vorticity  and 
stream- function  contours  are  presented  in  Fig.  C.3.  It  can  be  observed  that  the  grid  used  for  the 
calculation  is  quite  fine  and  should  provide  adequate  spatial  resolution.  The  grid  is  generated  using 
a  Laplacian  generation  system,  with  orthogonality  enforced  on  the  boundaries.  A  second  simulation 
was  performed  with  identical  parameters,  but  with  the  inlet  and  the  outlet  locations  placed  50  percent 
further  upstream  and  downstream  (-7.5  <  x  <  15.0  vs.  -5  <  x  <  10),  respectively,  and  the  grid 
increased  to  (675  x  45),  to  maintain  the  same  grid  spacing  as  for  the  shorter  configuration.  For 
subsonic  flows,  such  as  those  studied  in  the  present  work,  the  importance  of  the  placement  of 
inflow/outflow  boundaries  is  well  known.  However,  as  no  visual  difference  is  observed  in  the  results 


of  these  two  calculations,  it  can  be  concluded  that  the  solution  domain  of  the  shorter  configuration 
is  sufficiently  long.  Next,  simulation  results  were  obtained  for  the  case  of  a  (300  x  30)  initial  grid 
and  are  shown  in  Fig.  C.4.  Comparing  these  results  with  the  fmer-grid  results  shown  in  Fig  C.3, 
it  can  be  observed  that  the  separation  streamline  is  not  as  smooth  for  the  coarser  (300  x  30)  grid  as 
for  the  (450  x  45)  grid  case.  Further,  it  is  noted  that  the  reattachment  point  is  not  as  far 
downstream,  nor  is  it  as  sharply  defined  as  for  the  fine  grid  case.  Hence,  grid  adaption  was  then 
activated  and  the  simulation  continued  in  time.  Figure  C.5  clearly  shows  the  grid  adaption,  as  well 
as  the  smoothness  of  the  separating  streamline  and  the  regaining  of  the  separation  bubble  length. 

Two  simulations  with  fewer  grid-points  were  also  performed  in  order  to  test  the  procedure 
for  yet  coarser  grids.  Figures  C.6  and  C.8  present  results  for  (100  x  20)  and  (50  x  15)  fixed  grids, 
respectively.  The  results  were  difficult  to  obtain;  very  small  convergence  tolerances  were  required 
for  the  stream-function  equation  throughout  the  calculation.  Stalled  convergence,  with  oscillatory 
behavior  of  vorticity  residuals,  was  often  encountered,  necessitating  an  adjustment  of  the  time  step 
and  a  restart  of  the  calculation.  Then,  for  both  of  these  grids,  the  adaption  procedure  was  activated 
for  100  time  steps  and  the  simulation  advanced  to  steady  state.  The  results  obtained  using  the 
adaptive  technique  are  presented  in  Figs.  C.7  and  C.9  for  (100  x  20)  and  (50  x  15)  adapted  grids, 
respectively.  In  both  instances,  the  convergence  behavior  was  improved,  with  the  adapted  grid 
eliminating  the  need  for  user  intervention.  The  accuracy  of  the  simulation  improved  markedly  with 
the  adaption.  Separation  and  reattachment  points  became  more  sharply  defined  and  the  separation 
streamline  became  more  smooth.  The  vorticity  contours  show  no  visible  difference  for  these 
calculations,  and  hence,  are  not  repeated  for  each  grid  size.  The  axial  distribution  of  vorticity  is 
presented  for  two  grid  sizes,  namely,  (300  x  30)  and  (50  x  15)  in  Fig.  C.10.  The  variation  due  to 
grid  size  is  greatest  in  the  vicinity  of  the  separation  point.  The  length  of  the  separation  bubble 
gradually  decreases  as  the  number  of  grid  points  is  decreased,  reflecting  the  decreased  resolution. 
The  (100  x  20)  grid  size  provides  satisfactoiy  results;  the  contour  plots  presented  differ  only  slightly 
from  the  larger  (300  x  30)  adapted-grid  results,  in  the  region  immediately  surrounding  the 
reattachment  point.  The  (50  x  15)  grid  size  is  believed  to  provide  insufficient  resolution  even  with 
the  adaption,  as  the  separation  streamline  is  non-smooth.  However,  the  adapted  grid  is  believed  to 
provide  near  optimal  resolution  for  the  available  number  of  grid  points. 

The  same  geometric  configuration  was  also  studied  with  a  flow  Reynolds  number  of  1000. 
Simulations  for  two  grid  sizes,  namely,  (300  x  30)  and  (100  x  20)  were  performed  and  the  results 
are  presented  in  Figs.  C.tl  and  C.  12,  respectively.  The  two  sets  of  results  agree  well  with  each 
other,  demonstrating  that  the  adapted  (100  x  20)  grid  is  adequate  for  this  flow.  The  calculated 
velocity  profiles  at  various  streamwise  locations  are  compared  with  the  experimental  data  of  Ahmed 
and  Giddens  (1983)  in  Fig.  C.13.  The  agreement  is  close.  For  this  flow  situation,  the  flow  is 
indeed  axisymmetric.  It  can  be  observed  that  for  the  steady-state  flow  results  presented  in  Figs.  C.  1- 
C.  12,  the  reattachment  point  is  one  of  the  most  sensitive  quantities  with  respect  to  grid  size  and  grid- 
point  distribution.  Table  1  summarizes  the  reattachment  lengths  calculated  and  compares  the  results 
with  experimental  data  for  various  configurations.  The  sensitive  quantity  of  reattachment  length 
clearly  benefits  from  the  increased  resolution  capability  of  the  adaptive  procedure.  It  is  observed 
that  the  relationship  between  flow  Reynolds  number  and  reattachment  length  for  this  range  of 
Reynolds  number  is  linear,  as  theoretically  predicted  (Smith  [1977]). 

In  summary,  results  for  steady  flow  have  been  demonstrated  to  be  grid  independent.  Various 


geometrical  configurations  with  ACR  varying  from  0.2  to  0.89  and  XSL  —  1.0,  2.0  and  4.0  were 
studied;  none  maintained  attached  flow  for  Reynolds  number  greater  than  200.  The  grid-adaption 
procedure  has  also  been  shown  to  re-distribute  grid  points,  such  that  relatively  coarse  grids  provide 
results  comparable  to  those  with  much  finer  non-adapted  grids.  Typically,  the  calculations  required 
400-600  iterations  for  convergence  of  the  vorticity  to  residuals  less  than  l  x  10"*°.  In  addition,  the 
grid  adaption  is  activated  for  only  100  of  these  iterations.  Hence,  the  overhead  for  these  calculations 
is  less  than  10  percent.  Further,  the  numerical  simulation  has  nearly  reproduced  the  available 
measured  velocity  profiles. 

Unsteady  Flow  Results 

Simulation  results  are  presented  for  two  different  flow  Reynolds  numbers  for  a  geometric 
configuration  with  XSL  =  2  and  ACR  =  0.75  employing  a  (600  x  85)  grid.  These  configurations 
exhibit  inherently  unsteady  flow.  For  Re  =  2000,  Figs.  C.14a-g  present  the  grid,  vorticity  and 
instantaneous  stream- function  contours  for  various  instants  of  time  in  the  region  of  the  obstruction, 
while  Figs.  C.15a-g  present  the  same  data  corresponding  to  the  region  immediately  downstream  of 
this  location.  The  non-uniform  adapted  grid  is  obvious  in  the  sequence  of  these  figures.  The 
adaption  is  most  apparent  in  the  normal  direction,  with  the  largest  clustering  occurring  near  the  shear 
layer,  as  expected.  The  clustering  in  the  normal  direction  gradually  decreases  with  distance 
downstream,  as  the  fluid  mechanics  processes  distribute  the  vorticity  more  uniformly  in  the  normal 
direction.  The  clustering  near  the  wall  to  resolve  the  boundary  layer  is  quite  pronounced. 

Grid  clustering  in  the  streamwise  direction  is  less  obvious,  but  analysis  of  the  results 
simulated  using  various  values  for  the  weight  function  coefficients,  and  hence,  various  grid 
distributions,  has  indicated  that  the  distribution  attained  is  appropriate.  The  clustering  is  clearly 
visible  in  the  rapidly  accelerating  and  turning  region  near  the  throat.  Following  this  region,  the 
streamwise  distribution  of  the  grid  is  nearly  uniform,  up  to  approximately  x  =  5,  as  vortical 
structures  appear  throughout  this  region.  In  the  region  5  <  x  <  8,  the  streamwise  distribution  is 
less  constant  with  time,  as  pairing  of  vortical  structures  often  occurs  in  this  region.  The  distribution 
varies  from  nearly  uniform  when  several  structures  are  present,  to  a  more  clustered  distribution  when 
these  structures  have  combined  to  form  one  larger  structure.  For  x  >  8,  the  streamwise  distribution 
is  much  less  uniform,  as  discrete  vortical  structures  have  developed  with  regions  of  relatively  smooth 
vorticity  between  them.  The  difference  between  this  region  and  the  region  immediately  upstream 
is  that  the  former  region  always  contains  discrete  structures  and,  hence,  a  non-uniform  grid. 

Figure  C.  16  presents  the  flow  results  for  the  region  near  the  outlet.  It  is  observed  that,  for 
flow  at  the  relatively  low  Re  of  2000,  viscous  diffusion  has  smoothed  the  gradients  of  the  vortical 
structures  as  they  convect  downstream.  The  reduced  vortical  intensity  of  the  flow  is  reflected  in  the 
streamwise  sparsity  of  the  grid  in  the  downstream  region. 

Simulations  were  also  performed  for  Re  =  10,000.  However,  the  time  step  required  for 
temporal  resolution  was  such  that  racxing  the  flow  tor  any  meaningful  intervals  of  time  was 
computationally  very  expensive.  Figure  C.  17  shows  the  flow  results  for  the  region  near  the  outlet. 
For  Re  =  10,000,  the  flow  structures  convect  quickly  downstream,  and  hence,  maintain  their 
intensity  for  a  greater  streamwise  distance.  This  is  reflected  in  the  finer  spacing  of  the  grid  in  the 


downstream  region.  It  should  be  noted  that  the  results  presented  for  this  configuration  with  Re  = 
10,000  are  not  completely  resolved  in  time,  particularly  in  the  region  near  the  shear  layer.  In 
addition,  it  is  believed  that  the  axisymmetric  assumption  is  no  longer  appropriate  at  this  Reynolds 
number.  Currently,  work  is  being  undertaken  in  order  to  successfully  simulate  these  higher  Reynolds 
number  flows. 


2.C.7  Conclusion 

An  automated  flow-adaptive  grid  procedure  for  efficient  use  of  grid  points  to  resolve 
temporally  varying  regions  of  high  vorticity  and  high  vorticity  gradients  and  curvature  has  been 
developed.  The  technique  is  based  on  the  well  known  elliptic  grid-generation  procedures  and 
advances  the  grid  in  time  along  with  the  flow  solution,  as  the  current  grid  equations  are  parabolic 
in  nature.  Therefore,  the  present  technique  has  the  advantage  of  being  computationally  more 
efficient  than  a  purely  elliptic-equation-based  technique.  In  fact,  the  grid-transport  equations  fit  well 
into  unsteady  viscous  flow  solution  procedures.  The  procedure  uses  the  most  recent  flow  information 
to  adapt  the  grid,  not  information  from  previous  time  levels,  and  does  not  require  interpolation.  The 
relatively  small  overhead  associated  with  grid  migration  in  the  present  method  makes  the  technique 
computationally  practical  for  truly  unsteady  flow  problems.  The  relative  derivative  contribution  to 
the  weight  function  prevents  regions  of  large  vorticity  magnitude  from  swamping  regions  of  lower 
vorticity  which  may  also  require  resolution.  This  is  necessary  in  order  to  resolve  structures 
downstream  of  the  obstruction,  where  viscous  diffusion  has  diminished  the  gradients  appearing  in 
these  structures.  Finally,  many  of  the  existing  codes  employing  elliptic  grid-generation  systems  may 
be  modified  to  include  the  current  adaptive  grid  technique.  The  equations  combine  some  of  the  most 
desirable  features  of  several  types  of  grid  generators,  without  added  complexity  or  computational 
expense.  They  exhibit  the  efficiency  and  flow-tracking  capability  of  parabolic  equations,  while 
maintaining  the  smoothness  of  computationally  expensive  elliptic  equations. 

The  use  of  the  full  unsteady  equations  combined  with  the  moving  grid  capability  make  the 
code  very  versatile,  allowing  the  straightforward  implementation  of  procedures  to  handle  pulsatile 
flow  and/or  deforming  boundaries.  These  capabilities  are  essential  in  studies  of  flow-structure 
interactions  as  well  as  passive  and  active  flow  control.  The  capabilities  developed  here  are  also  of 
interest  to  the  biomechanics  community  for  the  study  of  flow  through  flexible  configurations. 
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JOURNAL  PAPERS 
PUBLISHED  AND  IN 
PREPARATION 


Books  and  Monographs 

Ghia,  K.N.,  and  Ghia,  U.,  "Computational  Fluid  Dynamics,"  section  in  1992  Yearbook  of 
Science  and  Technology.  Publisher  S.P.  Parker,  McGraw  Hill  Inc.,  NY. 

Celik,  I.,  Kobayashi,  T.,  Ghia,  K.N.,  Editors:  Advances  in  Numerical  Simulation  of 
Turbulent  Flows.  FED-Vol.  117,  ASME,  June  1991. 

Journal  Articles  and  Papers 
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Fluid  Dynamics  Research.  December,  1992. 
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Mechanics.  Editors:  R.C.  Baha  and  B.F.  Armaly,  Vol.  16,  October  1991,  pp.  28 1  - 
282. 
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of  Cincinnati,  Cincinnati,  <jhio,  January  1991. 
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Osswald,  G.A.,  Ghia,  K.N,  and  Ghia,  U.,  "Flow  Over  Arbitrary  Axisymmetric  Multi-Body 
Configurations  Using  Direct  Navier-Stokes  Simulations,"  Lecture  Notes  in  Physics, 
Vol.  371,  pp.  215-218,  September  1990. 

Osswald,  G.A.,  Ghia,  K.N.  and  Ghia,  U.,  "Direct  Simulation  of  Unsteady  Two-Dimensional 
and  Axisymmetric  Incompressible  Internal  Flows  Using  Navier-Stokes  Equations," 
Proceedings  of  IMACS  1st  International  Conference  on  Computational  Physics, 
University  of  Colorado  at  Boulder  Boulder,  Colorado,  pp.  252-262,  June  1990. 
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Ghia,  K.N.  and  Ghia,  U.,  "Simulation  Characterization  and  Control  of  Forced  Unsteady 
Viscous  Flow  Using  Navier-Stokes  Equations,"  Aerospace  Engineering  and 
Engineering  Mechanics  Report  AFL  91-11-76,  University  of  Cincinnati,  Cincinnati, 
Ohio,  November  1991. 

Ghia,  K.N.  and  Ghia,  U.,  "Direct-Solution  Techniques  for  Viscous  Flows  and  Their  Control," 
Department  of  Aerospace  Engineering  and  Engineering  Mechanics,  Report  AFL  91-5- 
75,  University  of  Cincinnati,  Cincinnati,  Ohio,  May  1991. 
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Section  5 


SCIENTIFIC 

INTERACTIONS-SEMINARS 
AND  PAPER  PRESENTATIONS 

Professional  Presentations: 

Invited  Lectures: 

Ghia,  K.N.,  Overview  Lecture  presented  at  QAI  Workshop  of  Fluid  Dynamics  anH 

Propulsive  Systems  Focus  Group  on  Needs  Assessment  Sandusky,  OH  March  6-7 
1992. 

Ghia,  K.N.,  presented  at  University  of  Tokyo.  Tokyo,  Japan,  August  30,  1991. 

Ghia,  U.,  presented  at  University  of  Tokyo  Tokyo,  Japan,  August  30,  1991. 

Ghia,  K.N.,  presented  at  Agronautical  .Development  Agency r  Bangalore,  India,  August  16 
1991. 

Ghia,  U.,  presented  at  Aeronautical  Development  Agency.  Bangalore,  India  August  16 
1991. 

Ghia,  K.N.,  presented  at  University  of  Taiwan.  Taipei,  Taiwan,  August  13,  1991. 

Ghia,  U.,  presented  at  .University  of  Taiwan.  Taipei,  Taiwan,  August  13,  1991. 

Ghia,  K.N.,  presented  at  National  Chenkung  University.  Tainan,  Taiwan,  August  12,  1991. 

Ghia,  U.,  presented  at  National  Chenkung  University.  Tainan,  Taiwan,  August  12,  1991. 

Ghia,  K.N.,  presented  at  ICASE/NASA  Langlev  Workshop  on  Transition  and  Turbulence 
Hampton,  V A  July  1991. 

Ghia,  U.  and  Osswald,  G.A.,  presented  at  General  Electric  Company.  Evendale  Ohio  March 
19,  1991. 

Ghia,  U  ,  Seminar  presented  at  Department  of  Aeronautical  and  Astronautical  Engineering 
Ohio  State  University,  Columbus,  Ohio,  March  6,  1991 .  . ^ 
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Ghia,  K.N,,  presented  at  NSF  Workshop  on  Advances  in  Computational  Methods  for 
Transport  Phenomena.  University  of  Kentucky,  Lexington,  Kentucky,  January  7-9, 
1991. 

Ghia,  K.N.,"  presented  at  The  Institute  of  Applied  Mathematics  and  Scientific  Computing. 
Indiana  University  -  Bloomington,  Bloomington,  Indiana,  November  1990. 

Invited  Papers; 

Ghia,  K.N.,  "Analysis  of  Two-Dimensional  Flow  Structure  of  Forced  Unsteady  Separated 
Flows,"  Leadoff  talk  presented  at  AFOSR  Workshop  on  Supermaneuverabilitv 
Physics  of  Unsteady  Flows  Past  Lifting  Surfaces  at  High  Angle  of  Attack."  Lehigh 
University,  Bethlehem,  PA,  April  9-10,  1992. 

Ghia,  U.,  ","  Leadoff  talk  presented  at  AFOSR  Workshop  on  Supermaneuverabilitv:  Physics 
of  Unsteady  Flows  Past  Lifting  Surfaces  at  High  Angle  of  Attack."  Lehigh  University, 
Bethlehem,  PA,  April  9-10,  1992. 

Ghia,  K.N.,  Yang,  J.,  Ghia,  U.,  Osswald,  G.A.,  "Analysis  of  Dynamic  Stall  Phenomenon 
Through  Simulation  of  Forced  Oscillatory  Flows,"  presented  at  Fifth  Symposium  on 
Numerical  and  Physical  Aspects  of  Aerodynamic  Flows.  California  State  University, 
Long  Beach,  CA,  January  13-15,  1992. 

Osswald,  G.A.,  Ghia,  K.N.  and  Ghia,  U.,  "Study  of  Forced  Unsteady  Separation  and  Its 
Control,"  presented  at  22nd  Midwestern  Mechanics  Conference.  University  of 
Missouri-Rolla,  MO,  October  6-9,  1991. 

Ghia,  K.N.  and  Rohling,  T.,  "Dynamical  Systems  Characterization  of  Unsteady  Flow,” 
presented  at  Fourth  International  Symposium  on  CFD.  University  of  California,  Davis, 
Davis,  CA,  September  1991. 

Ghia,  U.  and  Huang,  Y.,  "Evolution  and  Transport  of  Streamwise  Vorticity  in  Incompressible 
Flow,"  presented  at  Fourth  International  Symposium  on  CFD.  University  of  California, 
Davis,  Davis,  CA,  September  1991. 

Ghia,  K.N.,  Yang,  J.,  Osswald,  G.A.  and  Ghia  K.N.,  "Analysis  of  Structure  of  Force 
Unsteady  Separated  Flows  and  Its  Correlation  to  Dynamic  Stall  Vortex  Using 
Vorticity  Dynamics,"  presented  at  4th  International  Symposium  on  Supercomputing 
and  Experiments  in  Fluid  Dynamics.  Nobeyama,  Japan,  September  1991. 

Ghia,  U.,  Thornburg,  H.J.  and  Ghia,  K.N.,  "A  Flow-Adaptive  Grid  Technique  for  Simulation 
of  Time-Dependent  Vortical  Flow,"  presented  at  Fourth  Nobevama  Workshop  on 
Supercomputing  and  Experiments  in  Fluid  Dynamics.  Nobeyama,  Japan,  September  3- 
5,  1991. 

Ghia,  K.N.  and  Ghia,  U.,  "Three-Dimensional  Unsteady  Navier-Stokes  Analysis  and  Some 
Simplified  Analyses  for  Inlet  and  Nozzle  Flows,"  presented  at  OAI  Workshop  on  inlets 
and  Nozzles.  The  Ohio  Aerospace  Institute,  Brook  Park,  Ohio,  July  1990. 
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Ghia,  K.N.  and  Ghia,  U.,  "Direct  Numerical  Simulation  of  Dynamic  Stall  Phenomenon  Using 
Two-Dimensional  Unsteady  Navier-Stokes  Equations,”  IMACS  1st  Annual 
International  Conference  on  Computational  Physics.  Boulder,  Colorado,  June  1990 

Ghia,  K.N.,  "Chaotic  Flows  in  Open  System,"  IMACS  1st  Annual  International  Conference  on 
Computational  Physics.  Boulder,  Colorado,  June  1990. 

Ghia,  U.,  "Numerical  Experiments  on  Flow  Control,"  IMACS  1st  Annual  International 
Conference  on  Computational  Physics.  Boulder,  Colorado,  June  1990. 

Ghia,  U.  and  Ghia,  K.N.,  "Low  Speed  Forced  Unsteady  Flows,"  presented  at  International 
Symposium  on  Nonsteadv  Fluid  Dynamics.  ASME  Fluids  Engineering  Division  Meeting, 
Toronto,  Canada,  June  1990. 

Papers: 

Beltz,  S.,  Ghia,  K.N.  and  Rohling,  T.,  "Study  of  the  Effect  of  Camber  on  Window  of  Chaos  for 
Flow  Past  a  Joukowski  Airfoil,"  presented  at  AIAA  18th  Annual  Mini-Symposium  on 
Aerospace  Science  and  Technology.  Dayton,  Ohio,  March  1992. 

*Blodgett,  K.  and  Ghia,  K.N.,  "A  Spectral  Method  for  Solution  of  the  Navier-Stokes 
Equations  for  Arbitrary  Geometries,"  presented  at  AIAA  18th  Annual  Mini-Svmposium 
on  Aerospace  Science  and  Technology.  Dayton,  Ohio,  March  1992. 

Ding,  Z  ,  Ghia,  K.N.,  Osswald,  G  A.  and  Ghia,  U.,  "Simulation  of  Flow  Inside  a  Backstep 
Duct  Using  3-D  Unsteady  Incompressible  Navier-Stokes  Equations,"  presented  at 
AIAA  18th  Annual  Mini-Svmposium  on  Aerospace  Science  and  Technology.  Dayton. 
Ohio,  March  1992. 

Huang,  Y,,  Ghia,  U.  and  Osswald,  G.A.,  "3-D  Vorticity-Velocity  Simulation  of  an  Afterbody 
Flow,"  presented  at  AIAA  18th  Annual  Mini-Svmposium  on  Aerospace  Science  and 
Technology.  Dayton,  Ohio,  March  1992. 

Melde,  P.F.  and  Ghia,  K.N.,  "Investigation  of  Vorticed  Flows  Over  Slender  Delta  Wings," 
presented  at  AIAA  18th  Annual  Mini-Symposium  on  Aerospace  Science  and 
Technology,  Dayton,  Ohio,  March  1992. 

Noll,  C.,  Ghia,  K.N.  and  Muralirayan,  K  S  ,  "Analysis  of  Dynamic  Stall  Phenomenon  for 
Oscillating  Airfoil,"  presented  at  AIAA  18th  Annual  Mini-Svmposium  on  Aerospace 
Science  and  Technology.  Dayton,  Ohio,  March  1992. 

Wu,  H.C.,  Ghia,  K.N  and  Ghia,  U  ,  "Evaluation  of  Third-Order  Upwind  Difference  Schemes 
Using  Flow  Inside  Backstep  Channel,"  presented  at  AIAA  18th  Annual  Mini- 
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Symposium  on  Aerospace  Science  and  Technology.  Dayton,  Ohio,  March  1992 


Yang,  J.,  Ghia,  K.N.,  Osswald,  G.A.  and  Ghia,  U ,  "Active  Control  of  the  Dynamic  Stall 
Vortex  Using  Modulated  Suction  and  Injection,"  presented  at  AIAA  18th  Annual  Mini- 
Svmposium  on  Aerospace  Science  and  Technology,  Dayton,  Ohio,  March  1992 

Ghia,  K.N.,  Yang,  J.,  Osswald,  G.A.  and  Ghia,  U.,  "Study  of  the  Role  of  Unsteady  Separation 
in  the  Formation  of  Dynamic  Stall  Vortex,"  AIAA  Paper  92-0196,  AIAA  30th 
Aerospace  Sciences  Meeting.  Reno,  Nevada,  January  1991. 

Blodgett,  K.,  Ghia,  K.N.,  Osswald,  G.A.  and  Ghia,  U.,  "Pulsating  Flow  Past  an  Elliptic 
Cylinder,"  Bull.  Am.  Phvs.  Soc..  Vol.  36,  No.  10,  pp.  2656-2657,  1991,  presented  at  the 
44th  Annual  APS  Meeting.  Scottsdale,  Arizona,  November  1991. 

Ding,  Z.,  Ghia,  K.N.,  Osswald,  G.A.  and  Ghia,  U.,  "Simulation  of  Flow  Inside  a  Backstep 
Duct  by  Unsteady  Three-Dimensional  Navier-Stokes  Equations,"  Bull.  Am.  Phvs  Soc  . 
Vol.  36,  No.  10,  pp.  2621-2622,  1991,  presented  at  the  44th  Annual  APS  Meeting. 
Scottsdale,  Arizona,  November  1991. 

Thornburg,  H.J.,  Ghia,  U.  and  Ghia,  K.N.,  "Study  of  Free  Shear  Layer  Instability  and  Vortex 
Development  Using  a  Flow  Adaptive  Grid  Technique,"  Bull.  Am.  Phvs.  Soc..  Vol  36, 
No.  10,  pp.  2645-246,  1991,  presented  at  the  44th  Annual  APS  Meeting.  Scottsdale, 
Arizona,  November  1991. 

Ding,  Z.,  Ghia,  K.N.,  Ghia,  U.  and  Osswald,  G.A.,  "Hybrid  Method  for  Solution  of  Three- 
Dimensional  Unsteady  Incompressible  Navier-Stokes  Equations,"  presented  at  AIAA 
17th  Annual  Mini-Symposium  on  Aerospace  Science  and  Technology.  Dayton,  Ohio. 
March  1991. 

Ghia,  K.N.,  Yang,  J.,  Osswald,  G.A.  and  Ghia,  U,,  "Study  of  Dynamic  Stall  Mechanism  Using 
Simulation  of  Two-Dimensional  Unsteady  Navier-Stokes  Equations,"  presented  at 
AIAA  17th  Annual  Mini-Svmposium  on  Aerospace  Science  and  Technology.  Dayton, 
Ohio,  March  1991. 

Huang,  Y.,  Osswald,  G.A.  and  Ghia,  K.N.,  "Application  of  a  Multi-Grid  Method  in  Numerical 
Solution  of  Unsteady  3-D  Navier-Stokes  Equations,"  presented  at  AIAA  17th  Annual 
Mini-Svmposium  on  Aerospace  Science  and  Technology.  Dayton,  Ohio,  March  1991 

Rohling,  T ,  Ghia,  K.N.,  Osswald,  G.A.  and  Ghia,  U.,  "Nonlinear  Shear  Layer  Interaction  in  a 
Chaotic  External  Flow,"  presented  at  AIAA  17th  Annual  Mini-Svmposium  on 
Aerospace  Science  and  Technology.  Dayton,  Ohio,  March  1991. 

Shirooni,  S,  Ghia,  U  ,  Ghia,  K  N.  and  Osswald,  G  A  ,  "Vortex-Ring  Interactions  in  an 
Axisymmetric  Combustor  Flow  Configuration,"  presented  at  AIAA  17th  Annual  Mini- 
Svmposium  on  Aerospace  Science  and  Technology.  Dayton.  Ohio,  March  1991 


**Thomburgh,  H.J.,  Ghia,  U  and  Ghia,  K.N.,  "Efficient  Adaptive-Grid  Simulation  of 
Unsteady  Separated  Internal  Flow,"  presented  at  AJAA  17th  Annual  Mini-Svmposium 
on  Aerospace  Science  and  Technology.  Dayton,  Ohio,  March  1991. 

Yang,  J.,  Ghia,  K.N.,  Osswald,  G  A.  and  Ghia,  U.,  "Study  of  Initial  Acceleration  on  the 
Development  of  Dynamic  Stall  Vortex  for  an  Airfoil  Pitching  at  Constant  Rate," 
presented  at  AIAA  17th  Annual  Mini-Svmposium  on  Aerospace  Science  and 
Technology.  Dayton,  Ohio,  March  1991. 

Ghia,  K.N.,  Yang,  J.,  Osswald,  G.A.  and  Ghia,  U.,  "Study  of  Dynamic  Stall  Mechanism  Using 
Simulation  of  Two-Dimensional  Unsteady  Navier-Stokes  Equations,"  presented  at  29th 
Aerospace  Sciences  Meeting.  Reno,  Nevada,  January  7-10,  1991. 

Ghia,  U.,  Shirooni,  S.,  Ghia,  K.N.  and  Osswald,  G.A.,  "Examination  of  a  Vortex  Ring 
Interaction  Phenomenon  in  an  Axisymmetric  Flow,"  presented  at  29th  Aerospace 
Sciences  Meeting.  Reno,  Nevada,  January  7-10,  1991. 

Ghia,  K.N.,  Ghia,  U.  and  Osswald,  G.A.,  "Characterisation  of  Dynamic  Stall  Phenomenon 
Using  Two-Dimensional  Unsteady  Navier-Stokes  Equations,"  Proceedings  of 
NASA/AFOSR/ARO  Workshop  on  Physics  of  Forced  Unsteady  Separation.  Moffett 
Field,  California,  April  1990. 

Ghia,  K.N.,  "Analyses  of  Self-Excited  Oscillatory  Flows,"  presented  at  NASA/AFOSR/ARO 
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STUDENT  DEGREE  THESES 
AND  DISSERTATIONS 


M.S.  Degree  Theses 

Blodgett,  K.E.J.,  "Unsteady  Separated  Flow  Past  an  Elliptic  Cylinder  Using  the  Two- 
Dimensional  Incompressible  Navier-Stokes  Equations,"  M.S.  Thesis,  Department  of 
Aerospace  Engineering  and  Engineering  Mechanics,  University  of  Cincinnati, 
Cincinnati,  Ohio,  March  1990. 

Rohling,  T.,  "Simulation  of  Chaotic  Flows  Past  Airfoils  at  High  Incidence  Using  the  Unsteady 
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Ph.D.  Degree  Dissertations 

Thornburg,  H.J.,  "An  Adaptive-Grid  Technique  for  Simulation  of  Complex  Unsteady  Flows," 
Ph  D.  Dissertation.  Department  of  Mechanical,  Industrial  and  Nuclear  Engineering, 
University  of  Cincinnati,  Cincinnati,  Ohio,  August  1991. 
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Section  7 


TECHNICAL  APPLICATIONS 

Of  the  various  CFD  analyses  developed,  some  were  of  direct  use  to  the  technical 
community.  Although  to  our  knowledge,  none  of  these  analyses  were  used  in  the  de¬ 
velopment  of  any  specific  hardware,  they  are  being  used  in  preliminary  design  studies 
by  analysts  in  the  industry.  Some  of  these  analyses  are  also  being  used  by  other 
researchers  at  governmental  laboratories  to  improve  their  analyses.  The  following  is 
a  list  of  the  CFD  analyses  and  the  organizations  using  them. 


Analysis 

•  Three-Dimensional  Unsteady  Navier 
-Stokes  Analysis  Using  Iterative 
Technique  with  Multigrid  Acceleration 

•  Flow- Adaptive  Grid  Technique  For 
Simulation  of  Unsteady  Vortical  Flows 

•  Flow  Visualization  Technique 
for  Unsteady  Viscous  Flows 


Organization 

Copeland  Corporation 
Sidney,  Ohio 

McDonnel  Douglas  Research  Center 
St.  Louis,  Missouri 

The  Ohio  Supercomputer  Center 
Columbus,  Ohio 
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Appendix  A 
Figures  for  Section  2.A 
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Figure  A.l:  Inertial/body-fixed  coordinate  systems  and  boundary  conditions  for  flow 
past  an  arbitrary  maneuvering  body. 


h  H 

1/4  c  ' 

h - - - H 


Figure  A  2:  Depiction  of  forces  and  moments  on  the  maneuvering  airfoil. 


Figure  A. 3:  Representation  of  inviscid  flow  past  symmetric  N  AC  A  0015  airfoil  at 
angle  of  attar';  (with  circulation),  in  physical  and  various  transformed  planes. 


73 


Figure  A. 4:  Depiction  of  various  zones  for  streamwise  and  normal  1-D  analytical 
cubic  spline  clustering  transformations. 
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Figure  A. 5:  Grid  distribution  for  a  NACA  0015  airfoil  with  (444  x  101)  mesh  points. 


Constant  Pitch-up  Motion 
NACA  0015  Airfoil,  Re  =  45,000 
Computational  Plane,  C-Grid  (444,101) 

Vorticity  Distribution  (ADI  Solution) 

(a)  Time  =  1.50  a  =  20.54° 


Inc.  by  10  to  80  then  by  400 
Max. :  1061.42  Min. :  -3262.17 


Infinity  Zone 


(b)  Time  =  3.30  a  =  36.58 

Inc.  by  10  to  80  then  by  400 

Max. :  1511.13  Min. :  -3483.83 

Infinity  Zone 

Figure  A. 6:  Effect  of  grid  stretching  on  far- field  solution  -  vorticity  contours  fa) 
a  =  20.54°,  (b)  a  =  36.58°. 
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Figure  A.7:  Effect  of  grid  stretching  on  far-field  solution  -  disturbance  stream  function 
distribution,  (a)  a  =  20.54°,  (b)  a  =  36.58°. 
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•  Metric  Coefficient  Distributions 

1  2 

•  Computational  Plane  (4  ,%  ) 
•Grid  Size  ;  (444,101) 


Figure  A. 8:  Distribution  of  metric  coefficients,  for  round-off  error  study  (a)  Gil,  (b) 
G22/3. 
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•  Dirichlet-Poisson  Problem  for  BGE  Solution 

1  2 

•  Computational  Plane  (£,,$) 

•  Truncation  Error  in  Space  ;  0((A^  ),(A^2  )) 
(A£ 1  j2  =  5.1E-6,  (A^2)2=  1.0E-4 

•  Grid  Size  ;  (444,101) 


NORMALIZED  F<5.0  DISTRI8UTICN 


Max.  Absolute  Round-off  Error  -  2.91E-7 
at  (6.93E-1,  2.20E-1),  (3.07E-1,  2.20E-1) 


3  256-7 
t  636-7 
i  276-7 
9  066-8 
7.266-8 
5.456-8 


8 

7  316  9 
9  856-13 


ID 

Itt  >H 


o  so  or* 

? 


Max.  Absolute  Round-off  Error  =  3.46E-7 
at  (3.12E-1,  4.10E-1),  (6.88E-1,  4.10E-1) 


Figure  A. 9:  Estimation  of  accumulated  round-off  error  distribution  due  to  grid 
stretching  in  Block-Gaussian  Elimination  solution. 
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Grid  Size  ;  (330,75) 


(444,101) 


(544,121) 


pip  ■ 
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y  .1  , 

>.  ■  -- 

-  .  i-  .; 

-rS'iSg 
■  ?:-W& 

r  i'< 

■MiZguff 

(a)  Time  =  3.5 


Max.;  1083.99 
Min. 1085.80 


Max.;  1072.88 
Min.  ;-1071.91 


Max.;  1065.97 
Min.  -.-1063.73 


(b)  Time  =  4.5 


Max.;  1081.38 
Min. ;- 1083.20 


Max.;  1069.51 
Min. 107 1.30 


Max.;  1050.69 
Min.  -.-1075.02 


(c)  Time  =  5.5 


Max.;  1077.16 
Min.  ;-l08Z17 


Max.;  1090.82 
Min.  ;-1045.72 


Max.;  1075.99 
Min. 1050.52 


(d)  .Time  =  6.5 


Max.;  1081.13 
Min. ;- 1072.96 
Inc.  by  4  to  80  then  40 


Max.;  1085.65 
Min.  ;-1055.44 


Max.;  1044.52 
Min. ;- 1084.1. 


figure  A.  10;  Effect  of  three  different  grid  distributions  on  asymptotic  How  solution 
at  a  =  0°,  for  Re  =  45,000,  dr'1'  =  0.2.  (a-d)  vorticity  contours  near  TE,  (e)  time 
history  of  lift  coefficient. 


Grid  Size;  (330,75) 


(444,101) 


(544,121) 


Figure  A.  11:  Effect  of  three  different  grid  distributions  on  flow  past  a  NACA  0015 
airfoil  for  Re  =  45,000,  a+  —  0.2  -  contours  of  instantaneous  vorticity. 
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(a)  Times 1. 401,  cc=14.8l° 

Confident  of  Pressure  Distribution 


Wail  Vortidty  Distribution  on  Suction 


Surfaca 


(b)  Time=2.70l,  0=29  71° 

100  1 

1  Coaffident  of  Pressure  Distribution 

4Ku 


f  Wail  Vort'city  Distribution 


on  Suction  Surfaca 


lo  U  Comparison  of  Computad  ^  Exp«,m^ 


Lilt  Coariloants 


% 


,  ,  w  4Q 

AogJ*  erf  Attack  (Oog.) 


Constant  Pitch-Up  Motion 
NACA  0015  Airfoil 
Re=45,000 

Experiment  (Walker  et  al.,  1985, 

Re=47,500 ) 

Computed  (Grid  Size) 

^  030,75) 

(444,101) 

(544,121) 


figure  A.  12:  Effect  nf  tu  ~  "  ~~  - - 

aid°l1  for  Re  «  45,000,  d+Z  0.2*% c**  <?1Stri,butions  on  past  a  NACA  0015 
surface  at  a  =  14  gio  /u  r  j-  ^p  '  dist;ribution  and  wail  vnrf  •*.  ^ 

«  =  29.71°,  fr\  />'  d’.  .  p  ‘  dlstr‘bution  and  wail  vortidtv  *  7  °a  Suction 

l  -/  -a  -  distribution.  °rticity  on  suction  surface  at 


ONERA  Experiment 
WERLfe  (1976) 


Vorticity  Contours 


PRESENT 


Vorticity  Contours 


Velocity  Vectors 
MEHTA  (1977) 


Figure  A.  13:  Comparison  with,  numerical  results  of  Mehta  (1977)  and  experimental 
data  of  Werle  (1976)  for  NACA  0012  airfoil,  Re=5,00U,  k=l.O. 


PRESENT 


MEHTA  (1977) 


Figure  A.  14:  Comparison  with  numerical  results  of  Mehta  (1977)  for  NACA  0012 
airfoil,  Re=5,0QQ,  k=1.0  -  instantaneous  stream  function. 


PRESENT 


(a) 


(b)  Time  =  2.601,  a  =  18.57°  (Pitch-Up)  Time  =  4.601,  a  =  1 1.1 1°  (Pitch-Down) 


Figure  A.  16:  Comparison  with  numerical  results  of  Mehta  (1977)  and  Sankar  et  al. 
(1980)  for  NACA  0012  airfoil;  Re=5,000,  k=i.O.  (a)  vs.  a,  (b)  Cp  -  distribution 
for  a  =  18.57°  (pitch-up)  and  a  =  11.11°  (pitch-down). 
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Time  History  of  Lift  Coefficient  at  a  =  0° 


Instataneous  Vorticity  Contours  at  a  =  10.23° 
Start  to  pitch-up  at  (b)  w  .  637  « 


Min.  :  -2097.82 
CL  :  0.980 


Start  to  pitch-up  at  (c) 

Max.  :  S48.50 
Min.  :  -2083.26 
CL  :  0.990 


Start  to  Diteh-uo  at  ftf) 

- - “ — ^  Max. :  498.01 

Min.  :  -2166.54 

CL  :  1.106 


•  Constant  pitch-up  Motion 
NACA  0015  Airfoil 
Re  =  45,000 


Comparison  of  Lift  Coefficient 


Figure  A. 17:  Effect  of  three  different  initial  states  on  flow  past  a  NACA  0015  airfoil; 
Re  =  45,000,  a+  =  0.2  -  vorticity  contours,  lift  coefficients  and  unsteady  circulation. 
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a  =  24° 


a  =  27 


a  =  36° 


a  =  40° 


a  =  47° 


Figure  A. 18:  Comparison  with  experimental  data  of  Walker  et  al. ( 1 985)  for  flow  past 
a  NACA  0015  airfoil;  Re  =  45,000,  =  0.2  -  streaklines. 
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(a) 


4  0  Comparison  of  Computed  and  Experimental  Lift  Coefficients 


NACA  0015  Airfoil,  Re-45.000 


A 


Experiment  (Walker,  et  al.,  Re-47,500,  1985) 
Computed  (Visbal  and  Shang,  M-0.2,  1989) 
Computed  (Present,  M-0) 


20  30  40 

Angle  of  Attack  (Deg.) 


(b)  Time  =  1.401,  a  =14.81 


Time  =  2.701,  a  =  29.71 


Coefficient  of  Pressure  Distribution 
NACA  0015  Airfoil,  B— 15  000 


,,  .  - Emonmonl  (Waftor.  M  Bo-17  500,  tUS) 

-  - Computad 

9  V 


Figure  A.  19:  Comparison  with  experimental  data  of  Walker  et  al.(1985)  for  flow  past 
a  NACA  0015  airfoil;  Re  =  45,000,  q+  =  0.2.  (a)  CL  vs.  a,  (b)  Cp  -  distribution  for 
a  =  14.81°  and  a  =  29.71°. 
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Figure  A. 21:  Details  of  flow  structure  for  NACA  0015  airfoil;  Re=  45,000,  d+  =  0. 
-  instantaneous  vorticity  contours  for  a  =  22.83°, 26.27°, 28.56°,  and  29.71°. 
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Figure  A. 22:  Effect  of  suction/injection  for  NACA  0015  airfoil;  Re—  45,000,  and 
cC  =  0.2  -  unsteady  circulation,  velocity  vectors  and  Cp  -  distribution,  (i)  CASE  1: 
v,/Ux  =  0.045,  (ii)  CASE  2:  vJU„  =  0.035. 
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C/.\SE  1 


CASE  2 


Figure  A. 23:  Effect  of  suction/injection  for  NACA  0015  airfoil;  Re=  45,000,  and  a+  — 
0.2  -  instantaneous  vorticity  contours  and  wall  vorticity.  (i)  CASE  1:  v3/U00  =  0.045, 
(it)  CASE  2:  V'/Uoo  =  0.035. 


Comparison  of  Comguted  Lift  Coefficients 
NACA  0015  AUto«l.  fl#-4 5.000 


Wt ho*  Suci«ortftai«c»iafl 
CASE  I 
CASE  2 


20  0  _  240 

AnQis  at  A n*dt  (Dog.) 


Figure  A. 24:  Effect  of  suction/injection  for  NACA  0015  airfoil;  Re—  45,000,  and  a+  — 
0.2  -  Cl,Cd,Cm  and  Lq.  (i)  CASE  1:  u»/f/oo  =  0.045,  (ii)  CASE  2:  vs/U <*>  =  0.035. 


Velocity  Vectors 


Figure  A. 25:  Effect  of  suction/injection  for  NACA  0015  airfoil;  Re=  45,000,  or+  =  0.2 
and  va/Uoo  =  0.045  -  unsteady  circulation,  velocity  vectors  and  Cp  -  distribution,  (i) 
CASE  1:  t0  =  0.05,  (ii)  CASE  3:  t0  =  0.15. 
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CASE  1 


CASE  3 


Figure  A. 26:  Effect  of  suction/injection  for  NACA  0015  airfoil;  Re—  45,000,  q+  —  0.2 
and  v,/U0 0  =  0.045  -  instantaneous  vorticity  contours  and  wall  vorticity.  (i)  CASE 
1:  t0  =  0.05,  (ii)  CASE  3:  t0  =  0.15. 


(a) 


(b) 


(c) 


(d) 


Figure  A. 27:  Effect  of  suction/injection  for  NACA  0015  airfoil;  Re=  45,000,  a+  =  0.2 
and  v§/U„  =  0.045  -  CL,CD,CM  and  LD.  (i)  CASE  1:  t0  =  0.05,  (ii)  CASE  3:  ta  = 
0.15, 


101 


Appendix  B 

Tables  &  Figures  for  Section  2.B 
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\.  Or  id 

M«tnod 

25x25x25 

33x33x33 

43x40x^3 

MG-DGS 

0.4 

0.9 

2.5 

BGE 

32 

■i  ..  -  . 

135 

1000 

Table  B.l:  Storage  requirement  in  megawords  for  MG-DGS  and  BGE  methods  for 
velocity  problem 


U  -  Uniform  Grid  C  -  Clustered  Grid 


\ 

25x25x25 

33x33x33 

49x49x49 

U 

C 

u 

C 

u 

CPU* 

sac 

.39 

.58 

1.04 

1.31 

2.64 

wu 

12.8 

20.1 

16.2 

20.7 

13.5 

.41 

.57 

.48 

.58 

.42 

*  On  CRAY  Y-MP 


Table  B.2:  Some  results  for  model  problem  using  MG-DGS  method 


03 


^^Vlethod 

GricJ%!!!a^s=SiS5^ 

BGE 

MG-QGS0  ; 

■ 

25x25x25 

-4a 

0.5cxl0 

0.23x  10  ~4 

65x65x33 

-4  3  | 

0.22x10  ! 

s  Cyder  205  vectorized.  c  Uniform  Grid,  Pe*400,  T*  023,  t>  028.  £*i0 

0  CSAY-YMP  vectorized.  a  Clustered  Grid.  Re-3300.  T*18.  t*  01,  £-10 ' 


Table  B.3:  Comparison  of  CPU  time  (sec/per  grid/per  time  step)  for  BGE  and  MG- 
DGS  methods  for  cavity  flow  problem 


Table  4  Sizes  of  Secondary  Eddies  in  Symmetry  Plane 
at  T =15,  17  and  21. 


Eddy 

Time 

d3 

d. 

d, 

(Her  3  fle»30C0) 

T-15 

0 

.0652 

.3857 

T=17 

0 

.1063 

.3849 

.36' 

T-21 

,1082 

.2091 

.3831 

♦  Peed  from  Cn«rt 


Table  B.4:  Sizes  of  secondary  eddies  in  symmetry  plane  at  T  =  15,  17  and  21 
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mergence  history  for  DGS  scheme  for  model  problem  using  (25  x  25 
id. 


vergence  history  for  DGS  scheme  for  model  problem  using  (25  x  25 
y  clustered  grid. 
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Figure  B.3:  Convergence  history  for  DGS  scheme  for  model  problem  using  (25  x  25 
x  25)  strongly  clustered  grid. 


Figure  B.4:  Convergence  history  for  DGS  scheme  for  velocity  solution  of  shear-driven 
3-D  cavity  flow. 
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X 


Symmetry  Plane  MM 


Figure  B.5:  Configuration  of  3-D  cavity  flow. 


Figure  B.6:  Flow  pattern  in  symmetry  plane  MM. 


W  (M 

Figure  B.7:  Visualization  of  Ref.  [B.14]  for  cavity  flow  (spanwise  aspect  ratio  3;1); 
(a)  symmetry  plane  MM;  (b)  lateral  plane  NN  (i1  =  0.8).  Re  =  3200,  time-averaged. 
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(•)  (>>} 

Figure  B.8:  Simulation  results  for  cavity  flow  (spanwise  aspect  ratio  3:1)  in  symmetry 
plane  MM;  (a)  velocity  vectors;  (b)  vorticity  contours.  Re  =  3300,  T  =  15,  17,  21. 


Figure  B.9:  Simulation  results  for  cavity  flow  (spanwise  aspect  ratio  3:1)  in  lateral 
plane  NN;  (a)  velocity  vectors;  (b)  vorticity  contours.  Re  =  3300,  T  =  15,  17,  21. 
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Figure  B.10:  Tangential  velocity  vector  plots  and  normal  vorticity  contour  in  sym¬ 
metry  plane  z  =  0  for  t  =  20,  25,  30,  35,  45,  50. 
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Figure  B.12:  No 
and  x  =  0.265  fo 


Appendix  C 

Table  &  Figures  for  Section  2.C 
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|  Source  Rea ttachroent  point  location  Xr  S 

|  ACR  »  0.50.  XSL  -  2.0  jj 

1  RE  •  500  | 

Exp«flm«nU)  daU  of  Ahmad  and  Giddana  (1983) 

■  4.5 

Preaant  Computation,  (300  x  30)  Fixed  grid 

42 

Praaent  Computation,  (300  x  30)  Adapted  grid 

4.45 

Preaant  Computation,  (tOO  x  20)  Adapted  grid 

4.2 

Praaent  Computation,  (50  x  15)  Adapted  grid 

4.0 

i  n*  1 

. 

Experimental  data  of  Ahmed  and  G  Id  den  a  (1953) 

-55 

Preaant  Computation,  (300  x  30)  Adapted 

54 

|  R«  *_JQW  I 

Experimental  data  of  Ahmad  and  Glddene  (1983) 

8  0*  (Not  wall  defined) 

Praaent  Computation,  (300  x  30)  Adapted  grid 

89 

Praaent  Computation,  (100  x  20)  Adapted  grid 

6.7 

Table  C-l:  Comparison  of  reattachment  length.  ACR  =  0.50,  XSL  =  2.0 


t 


rv3 

e 


Figure  C.l:  Generalized  axisymmetric  coordinate  transformation  T. 
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Figure  C.2:  Schematic  of  geometric  configuration. 


R#  -  500  ,  (450x45)  Grid  ,  Tim*  -  32.333 


30  2.1  12  -0.3  0«  IS  24  33  42  5.1 

X 


3  0  2.1  -12  03  O.a  IS  2-4  33  4  2  5.1  5  0 

X 


INSTANTANEOUS  STREAM  FUNCTION  CONTOURS:  -  0  2523 


3  0  2.1  -1  2  -03  O.t  IS  2.4  33  4  2  Si 

X 


Figure  C.3:  Steady  flow,  50  %  constriction  by  area,  XSL  =  2.0,  (450  x  45)  fixed  grid, 
Re  =  500. 
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Re  *  500  .  (100x20)  Grid  ,  Time  -  50.0 


Figure  C.6:  Steady  flow,  50  %  constriction  by  area,  XSL  =  2.0,  (100  x  20)  fixed  grid, 
Re  =  500. 


Re  -  500  ,  (100x20)  Grid  .  Time  *  87.0 


Figure  C.7:  Steady  flow,  50  %  constriction  by  area,  XSL  =  2.0,  (100  x  20)  adapted 
grid,  Re  =  500. 
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Figure  C.8:  Steady  flow,  50  %  constriction  by  area,  XSL  =  2.0,  (50  x  15)  fixed  grid, 
Re  =  500. 
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Figure  C.10:  Axial  variation  of  wall  vortidty,  steady  flow,  50  %  constriction  by  area. 
Re  =  500. 
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Figure  C.ll:  Steady  flow,  50  %  constriction  by  area,  XSL  =  2.0,  (300  x  40)  adapted 
grid,  Re  =  1000. 


Re  *  1000  .  (100x20)  Grid  ,  Time  -  139.0 


Figure  C.12:  Steady  flow,  50  %  constriction  by  area,  XSL  =  2.0,  (100  x  20)  adapted 
grid,  Re  =  1000. 


R/R 


IMJK 


R/R 


udflyc 


Figure  C.13:  Axial  velocity  profiles,  50  %  constriction  by  area,  XSL  =  2.0,  (300  x 
30)  adapted  grid,  Re  =  1000. 
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Re  =  2000  ,  (600x85)  Grid  .  Time  -  50.2 
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b.T  -  50.2. 

Figure  C.14:  Unsteady  flow  adaption,  throat  region,  Re  =  2000,  XSL  =  2.0,  ACR 
0.75,  (600  x  85)  grid. 
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Re  »  2000  ,  (600x85)  Grid  ,  Time  -  50.4 
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Figure  C.14:  Continued. 
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Re  -  2000  .  (600x85)  Grid  ,  Time  *  51.0 
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Figure  C.14:  Continued. 
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g.  t  *  51.2. 

Figure  C.14:  Continued. 


Re  -  2000  .  (600x85)  Grid  ,  Time  -  50.0 
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a.  x  »  50.0. 

Figure  C.15:  Unsteady  flow  adaption,  downstream  region,  Re  =  2000,  XSL  =  2.0, 
ACR  =  0.75,  (600  x  85)  grid. 
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Re  -  2000  .  (600x85)  Grid  ,  Time  *  50.2 
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Re  -  2000  ,  (600x85)  Grid  ,  Time  =  50.4 
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Figure  C.15:  Continued. 
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R*  -  2000  .  (600x06)  Grid  ,  Tin*  »  50.6 


d.  t  -  50.6. 


Re  *  2000  .  (600x85)  Grid  .  Time  -  50.8 


e .  t  *  50.8. 


;  Figure  C  l  5:  Continued. 
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Figure  C.  15:  Continued. 
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